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Abstract 

The foundation of the local energy-density functional method to describe the nuclear 
ground-state properties is given. The method is used to investigate differential ob- 
servables such as the odd-even mass differences and odd-even effects in charge radii. 
For a few isotope chains of spherical nuclei, the calculations are performed with an 
exact treatment of the Gor'kov equations in the coordinate-space representation. 
A zero-range cutoff density-dependent pairing interaction with a density-gradient 
term is used. The evolution of charge radii and nucleon separation energies is re- 
produced reasonably well including kinks at magic neutron numbers and sizes of 
^ ■ staggering. It is shown that the density-dependent pairing may also induce sizeable 

staggering and kinks in the evolution of the mean energies of multipole excitations. 
The results are compared with the conventional mean field Skyrme-HFB and rela- 
tivistic Hartree-BCS calculations. With the formulated approach, an extrapolation 
from the pairing properties of finite nuclei to pairing in infinite matter is considered, 
and the dilute limit near the critical point, at which the regime changes from weak 
to strong pairing, is discussed. 
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1 Introduction 



For a long time, the fine structure of the isotopic dependence, i.e. dependence 
on the neutron number, of nuclear charge radii could not be explained satis- 
factorily. Because of the apparent lack of relevant experimental information, 
the effective particle-particle (pp) force which leads to the observed pairing 
properties of nuclei, has mostly been assumed in a very simple form, just suf- 
ficient to produce reasonable gap parameters A extracted from the observed 
odd-even staggering in the nucleon separation energies. Even more sophisti- 
cated effective interactions, e.g. the Gogny force [1-3], do not include a density 
dependence in the pp channel. 

It then has been demonstrated in Hartree-Fock-Bogolyubov (HFB) type cal- 
culations that three- or four-body forces are essential to reproduce the exper- 
imental data on isotope shifts of nuclear charge radii [4-6], and indeed that 
these isotope shifts give indirect experimental information on the effective in- 
medium many-body force, or what is equivalent, on the density dependence 
of the effective interaction, particularly in the pp channel. 

It should be emphasized that a better knowledge of the density dependence 
of the nuclear pairing force would give the possibility to predict the pairing 
gap as a function of nuclear matter density which is, in particular, of great 
importance for understanding the pairing phenomena in neutron stars. It is 
well known that at present the pairing gap can not be obtained with sufficient 
accuracy from nuclear matter calculations based on bare NN interaction. The 
empirical information gained from the studies of real nuclei, specifically from 
the combined analyses of the nucleon separation energies and isotope shifts in 
charge radii, seems to be indispensable in this respect. A more general remark 
is that, as one may notice, a "good" microscopic theory, which could supply an 
effective interaction for describing the nuclear ground states and low-energy 
nuclear structure on a satisfactory level, is still lacking. The presently most 
successful simultaneous description of the bulk nuclear properties, such as 
binding energies and radii, throughout the periodic chart is achieved with 
phenomenological density- dependent interaction such as the Skyrme force [7]. 
This can be traced to the fact that some major effects produced by the den- 
sity dependence of the effective interaction in the particle-hole (ph) channel, 
which is the second variational derivative of the corresponding energy-density 
functional with respect to normal density, are now established fairly welQ On 
the same ground, one expects that simultaneous description of the differen- 
tial observables, such as odd-even mass differences and the odd-even effects in 



3 One may still notice that a universal and unique parametrization of the ph force 
has not yet been found. In particular, considerable effort is still continuing to opti- 
mize the HF part of the Skyrme-type functional [8,9]. 
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radii, would shed light on the density dependence of the effective interaction 
in the pp channel. 

In [4-6] an effective interaction has been chosen which consisted of two- and 
three- (or four-) body parts, with parameters adjusted independently of those 
of the single-particle potential used for the shell model description of the 
reference nucleus. Only differences with respect of this reference nucleus have 
been calculated in a self-consistent way. 

This procedure corresponds to the philosophy of the Landau-Migdal theory 
of finite Fermi systems [10], where there is no simple connection between the 
single particle well parameters and the effective interaction of quasiparticles 
close to the Fermi surface. However, there are consistency requirements [11-16] 
and since the parameters have been fitted to a large number of data, a quite 
reliable set of force parameters is available [16]. The interaction of Refs. [5,6] 
has been restricted by the requirement that it should reduce to the Migdal 
force in the ph channel, the three-body part leading to the density dependence. 

In the meantime, since Migdal had formulated his theory of finite Fermi sys- 
tems (FFS), there has been much progress in Hartree-Fock calculations [17- 
19,3], as well as in the self-consistent version of the FFS theory [15,20,16], 
demonstrating that initially Migdal may have been too pessimistic concerning 
the possibility to use the same interaction for ground state properties and for 
low-lying excited states. 

In deriving the HF or HFB equations, the starting point is the minimization 
of the expectation value of the energy, expressed by an effective interaction. 
The latter is considered as a substitute for the G-matrix derived from the true 
NN interaction, appearing in some kind of "effective Hamiltonian" . Therefore, 
ideally, the same interaction should be used in the pp (hh) and the particle- 
hole (ph) channel. The energy is then obtained from the interaction and the 
state vector which is assumed to be a Slater determinant. 

The self-consistent FFS theory or energy density functional (EDF) method 
starts from relations between the total energy, expressed as a functional of the 
density, the single particle potential, the effective interaction, and the quasi- 
particle density: Instead of parameterizing the interaction and deriving the 
energy from it in a certain approximation, an ansatz for the energy functional 
is assumed from which the other quantities can be derived. The interaction 
in the pp (hh) channel is obtained from the energy functional by a different 
procedure than that in the ph channel, and therefore these two interactions 
are different, which is of great practical importance. This feature is shared 
with Migdal's theory of finite Fermi systems. 

Otherwise, from a technical point of view, there is little difference between 
the EDF method and the HF or, with pairing, HFB method. The latter cor- 
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respond to special choices of the energy functional. On the other hand in the 
EDF method, to compute the density, a Slater determinant or HFB-type state 
vector is employed. Thus the differences which exist in principle, due to the 
necessary approximations are of little practical importance. 

In a recent paper [21] the EDF method has been applied to magic nuclei 40 ' 48 Ca 
and 208 Pb, investigating especially the influence of the spin-orbit interaction 
on collective states. Besides demonstrating the importance of including the 
two-body spin orbit force for the consistency of the method in describing the 
excited states, with the chosen set of parameters good results on the ground 
state properties have been obtained. 

Here an analogous EDF approach is applied to open-shell nuclei with special 
attention to the pairing part of the functional. First results obtained mainly in 
diagonal-pairing approximation (i.e. with a HF-BCS formalism) have already 
been given in a few short publications [22-24] (the paper [23] gives the first 
application of the EDF approach to deformed nuclei and contains the results 
of a HFB-type calculation for dysprosium isotopes). In the present paper we 
give the results of more elaborate calculations based on the coordinate-space 
technique developed in [25] which allows us, for a local pairing field A(r), to 
solve the Gor'kov equations in spherical finite systems without any approxi- 
mations [26]. The latter approach, even in the case of contact pairing force, 
corresponds, to a good approximation as will be shown in the present paper, to 
a full treatment of the HFB problem by applying the general variational prin- 
ciple to the EDF with a fixed energy cutoff e c > e F (e F is the Fermi energy). 
It is known from nuclear matter calculations that s-wave pairing vanishes at 
or slightly above the equilibrium density ~ 0.16 fm -3 ; and so, in our local 
EDF approach, the cutoff is chosen to be larger than the corresponding Fermi 
energy e F ~ 37 MeV. In actual calculations we shall use e c = 40 MeV. 

In Section 2 we give some theoretical foundation of the local EDF method for 
superfluid nuclei. In Section 3 the normal part of the energy-density functional 
with the parameters used in the calculations is described. In Section 4 the 
pairing part of the functional is presented. There, a possible extrapolation 
to uniform nuclear matter and the dilute limit with weak and strong pairing 
is also considered. Section 5 is devoted to the analysis of specific numerical 
results for a few isotopic chains of spherical nuclei and to the comparison of 
these results with experimental observations and other theoretical approaches. 
The influence of the density-dependent pairing on the evolution of the mean 
energies of multipole excitations is also studied using the self-consistent sum 
rule approach. In Section 6, the contribution of the ground state ("phonon") 
correlations to the charge radii is analyzed. In Section 7 the conclusions are 
summarized. Some important theoretical and technical aspects are presented 
in the Appendices A, B, and C. Appendix A contains a thorough formulation of 
the generalized variational principle for the systems with pairing correlations 
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which could be described by a local cutoff EDF. In Appendix B, the expression 
for the pairing energy is derived by using the Green's function formalism. In 
Appendix C, we give a detailed description of the coordinate-space technique 
which we apply here to spherical nuclei. 



2 Method of the local energy functional with pairing 

In this section we give some foundation of the local-energy-density functional 
method which we apply here. Let us first briefly mention a few approaches 
used to calculate the nuclear ground state energy which is the expectation 
value of the Hamiltonian H (with a "bare" NN-interaction) over the exact 
ground state vector |<3>): 

E=($\H\$). (2.1) 

In general, the vector |$) is a functional of the particle density g p [27] given 
by 

e P(f 1 ) = Tr rt <$|*t(i)«r(i)|$) j (2.2) 

where s and t are the spin and isospin indices and ^(1) and are the 
particle creation and annihilation operators, respectively, (1) = {{r, s z ,r 3 }i). 
Therefore, the energy is a functional of g p : 

E = E[g p ] . 

It is a difficult problem to construct such a functional and calculate the energy 
for many body systems. Formally, one can introduce single-particle degrees of 
freedom and, assuming that there is no pairing condensate, extract the single- 
particle motion on the HF level by writing^ 

E\p] = (13F\H\HF) + E caa \p], (2.3) 

where |HF) is the HF vacuum — a Slater determinant of single-particle wave 
functions fa. The latter are defined in a self-consistent manner by the equation 

Of course, in the nuclear case, this expression with bare H is meaningless since 
(HF|if|HF) becomes infinite due to the repulsive core. It is understood then that, 
for the first term in (2.3), the H is taken with some kind of microscopically derived 
effective interaction, for example with the Brueckner G-matrix (see, e.g., the dis- 
cussion in [28]). Because of inevitable approximations, such a replacement would 
never lead to the vanishing of the second term, E corv . 
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h(pi = in which the single-particle Hamiltonian h is given by the variational 
derivative of the energy functional with respect to the single-particle density 
matrix p: h = 5E/5p. The second term, E corr , is the dynamical correlation 
energy (short range correlations, RPA correlations, many-particle correlations, 
parquet diagrams, etc.). The density matrix p may be associated with the 
Landau-Migdal quasiparticle density matrix defined as 

p(l,2) = (HF|^(2)^(l)|HF>, 

with ifv (if)) the quasiparticle creation (annihilation) operators. It is in turn a 
functional of the particle density g p : 

p = p[Q»] . 

In homogeneous infinite nuclear matter, because of equality of the particle and 
quasiparticle numbers, the following relation between the two densities holds: 

£ p = p = Tr st p(l,l). (2.4) 

In an inhomogeneous system, the local difference between the densities could 
be attributed to the quasiparticle form factor [15]. Using the effective radius 
approximation one may write 

gP(f) = (1 + Ji^VXrO , (2.5) 
where R q is the effective quasiparticle radius. 

Hence the mean square nuclear radius for the particles may differ from that 
for the quasiparticles: 

(ry = {r 2 ) + R 2 q . (2.6) 

Some additional effects may be related with a possible change of the internal 
structure of nucleons in nuclear medium. Such effects are not included explic- 
itly in the ordinary approaches, based either on a "bare" Hamiltonian H or on 
some kind of effective interaction, in which the nucleons or quasiparticles are 
considered as point-like objects^. Within the EDF approach, however, tracing 
back to the Hohenberg-Kohn density-functional theory [27], if the functional 
were known, its minimization would determine both the ground-state energy 

5 To get the charge density distribution, and the charge radius of a nucleus, the 
point nucleon densities are usually folded with the free nucleon charge form factors. 
These may also change in the nuclear medium. 
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and the correct particle density g p (r). Anyway, for the differences between 
mean square radii of nuclei, which is of our main interest here, one has 



<5( r 2 )P = 5(r 2 ) . 



(2.7) 



The common problem of any many-body theory is how to calculate the dy- 
namical correlation energy E corr or how to take it into account in some effective 
way. In nuclear physics a few self-consistent methods are very popular: 

- The HF method with effective density- dependent interaction (e.g., with 
zero-range Skyrme [19] or finite-range Gogny force [3]). In this method 
the energy is calculated as a ground state expectation value of an effec- 
tive Hamiltonian H e g taken over a Slater determinant, 



- The relativistic mean field model [29,30] based on a phenomenological rel- 
ativists Lagrangian which includes mesonic and nucleonic degrees of free- 
dom. In most applications to the ground states of finite nuclei, the relativis- 
tic Hartree formalism with static meson fields, and with no-sea approxima- 
tion for nucleons, has been used (see [31] and references therein). 

- The quasiparticle Lagrangian method (QLM) which is an extension [15] of 
the Landau-Migdal quasiparticle concept. In this method an effective local 
quasiparticle Lagrangian is constructed taking into account the first-order 
energy-dependence of the nucleon mass operator and the requirements im- 
posed by self-consistency relations [13]. The QLM can be reformulated in 
terms of an effective Hamiltonian [32] in which the (linear) energy depen- 
dence of the nucleon mass operator is completely hidden so that this ap- 
proach becomes equivalent to the EDF method (see also the discussion of 
this point in [20]). 

- The method of the effective quasiparticle local energy density functional 
(which is referred to as EDF here). The possibility to use this method is 
based on the existence theorem of Hohenberg and Kohn [27]. With the 
Kohn-Sham quasiparticle formalism [33], this theorem allows one to write 
the nuclear EDF in the form 



where the first term is taken with the free kinetic energy operator t = p 2 /2m 
(m is the bare nucleon mass). The ground state energy is then determined 
by making this EDF stationary with respect to infinitesimal variations of the 
single-particle wave functions belonging to the class of the Slater determi- 
nants from which the density matrix p (and the density p) can be calculated 
self-consistently. This method is flexible in the sense that the EDF does not 



E[p] = (HF|# cff |HF>. 



(2.8) 



E\p]=Tr(tp) + EU[p], 



(2.9) 
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(and should not) expose the same symmetry properties for the effective 
force, which is the second variational derivative of the EDF with respect to 
p, as the underlying bare NN interaction or effective Skyrme-type forces. 
Particularly, the force in the ph channel obtained from the EDF does not 
come out to be antisymmetrized, and the effective force in the pp channel 
needed for describing the nuclear pairing properties has to be derived in 
a different way. As a result, the matrix elements of the pp force, although 
they should be antisymmetrized, have no direct connection, say by means of 
simple angular momentum recoupling like the Pandya transformation, with 
those of the ph forceQ This is, as discussed in the Introduction, in accord 
with the philosophy of Migdal theory of finite Fermi systems [10]. (For an 
example of the phenomenological nuclear EDF see [20,21]). 

In most of the applications of the conventional functionals to finite nuclei the 
pairing correlations are introduced either at the BCS or the HFB level (see, 
e.g., Ref. [37] and references therein). When the contact pairing interaction is 
used, the pairing part of the functional is usually evaluated within a truncated 
space of the quasiparticle levels imposing an energy cutoff which is often chosen 
with some freedom but not too far from the Fermi surface. A more rigorous 
approach should be based on the general variational principle for the cutoff 
functional which uses the same quasiparticle basis both in the ph and pp 
channel. 

We proceed as follows. In nuclei with pairing correlations one works with ap- 
proximate state vectors which are not eigenstates of the particle number. One 
assumes nonzero anomalous expectation values, 0(1,2) = ($1^(1)^1/(2) |$) ^ 
0, and the energy of a superfluid nucleus may be given by a functional of 
the generalized density matrix R containing both normal, p, and anomalous, 
0, components (see Appendix A). In analogy to eq. (2.3), the energy of a 



6 This was explicitly demonstrated, already on the level of the "ladder" approxima- 
tion, in [34]. We note the term "force", used to name the second variational deriva- 
tive of the EDF with respect to the density matrix, may be somewhat confusing in 
this context. In fact, this derivative, for infinite systems, has strict correspondence 
to the quasiparticle interaction amplitude on the Fermi surface introduced in the 
Landau theory of Fermi liquids [35] . The total quasiparticle scattering amplitude 
is, of course, antisymmetric due to the Pauli principle resulting in sum rules for 
the Landau parameters. These parameters should satisfy the Pomeranchuk stabil- 
ity conditions [36] in order to prevent the collapse of the system with respect to 
the low-energy ph excitations and provide the ground state energy to be a mini- 
mum rather than simply stationary. For finite systems, this means that all the RPA 
solutions obtained with the self-consistent ph interaction, taken as the second vari- 
ational derivative of EDF at the stationary point, should have real frequencies u> 
such that uj 2 > 0. 
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superfluid nucleus can be written as 

E[R] = (HFB|#|HFB) + E corr [R] (2.10) 

with |HFB) the HFB quasiparticle vacuum and E corr the dynamical correlation 
energy. 

The generalized density matrix R, in turn, is a functional of g p and therefore 

E = E[R[g p }] = E[g p ] . (2.11) 

The method we follow here is based on an effective quasiparticle energy func- 
tional of the generalized density matrix: 

E[R}= E kiQ [p} + E int [p,0], (2.12) 

where E khl [p] = Tr(ip) , and E int [p,i>] = E illt(noimal) [p] + £ a nomai[p, v\ . The 
anomalous energy -Eanomai is chosen such that it vanishes in the limit v — > 0. 

It should be emphasized that the weak-pairing approximation |A| <C ep is 
assumed throughout this paper. That means we need to retain only the first- 
order term ~ v 2 in the ^-dependent, anomalous part of the energy functional. 
One may then write 

E ml [p,u] = l(u^ pp [p]u) , (2.13) 

where P pp is an antisymmetrized effective interaction in the pp channel and 
the round brackets (...) imply integration and summation over all variables. 

To calculate the ground state properties, one can now use the general varia- 
tional principle with two constraints, 

(HFB | N(p) | HFB) = N(p) = N , (2.14) 

R 2 = R, (2.15) 

leading to the variational functional of the form 

I[R] = E[R] - pN(p) - TrA(R - R 2 ) , (2.16) 

where iV is the particle number, p the chemical potential, and A the matrix 
of Lagrange parameters (see, e.g. Ref. [28], and Appendix A for more details). 
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The major difficulty in implementing this general variational principle in prac- 
tice is connected with the anomalous part of the energy functional. The anoma- 
lous energy (2.13) can be calculated if one knows a solution of the gap equation 
for the pairing field A and the anomalous density matrix v. In general, the 
gap equation, 

A = 1^0, (2.17) 

is nonlocal and its solution (starting, for example, from a realistic bare NN 
interaction [38-41]), even for uniform nuclear matter, poses serious problems. 
First attempts were made very recently to construct an effective pairing in- 
teraction for semi-infinite nuclear matter [42] and for finite nuclei [43] with 
an approximate version of the Brueckner approach. Those studies are in an 
initial stage, and so far they do not give any guidance how to choose an ef- 
fective pairing interaction, in particular its density dependence, that could be 
used in nuclear structure calculations. It is assumed that a simple universal 
effective interaction in the pp-channel can be invented to correctly describe 
the nuclear pairing properties. Our intention is to show that, to a good ap- 
proximation in the case of weak pairing |A| <C eF, the EDF method and the 
general variational principle can be used with an effective density-dependent 
contact pp-interaction. 

We proceed with the formal development by using the Green's function for- 
malism as described in detail in Appendices A and B. Here we give only a 
brief account of the main issues. 

We introduce an arbitrary cutoff e c in the energy space, but such that e c > e-p, 
and split the generalized density matrix into two parts, 

R = R C + S C R, (2.18) 

where 5 C R is related to the integration over energies |e| > e c . 
The gap equation is renormalized to yield 

A = (2-19) 

where z> c is the cutoff anomalous density matrix and JFJ is the effective anti- 
symmetrized pp-interaction in which the contribution coming from the energy 
region |e| > e c far from the Fermi surface is included by renormalization (see 
Appendix A). 

For homogeneous infinite matter with weak pairing (|A| <C ep), it is shown 
that the variational functional E — fiN does not change in first order in |A| 2 /e F 
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upon variation with respect to 5 C R. The total energy of the system and the 
chemical potential also remain the same in first order in | A| 2 /ep if one imposes 
the particle number constraint (2.14) for the cutoff functional. To a good 
approximation, as discussed in Appendix A, this should also be valid for finite 
(heavy) nuclei. 

Such an outcome may be understood by noting that the major pairing effects, 
if |A| <C ep which is generally the case in nuclear matter and in finite nuclei, 
are developed near the Fermi surface. In this connection we mention the well- 
known fact that the pairing energy, in the BCS approximation, is defined by 
a sum concentrated near the Fermi surface (see, for example, Ref. [28]): 

£pair~-inA 2 , (2.20) 



with n the average level density and A the average energy gap in the vicinity 
of the Fermi surface. In infinite matter, this corresponds to the pairing energy 
per particle (see Appendix B): 

^ = _3A !M 

N 8 e F K J 



It follows that, with the cutoff functional, this leading pairing contribution to 
the energy of the system is exactly accounted for. 

Thus we find that nuclear ground state properties can be described by applying 
the general variational principle to minimize the cutoff functional which has 
exactly the same form as in eq. (2.16) with the constraint (2.14) but with R 
replaced by R c : 

I C [R] = E kin [p c ] + E°Jp c , v c ] - p c N c (p c ) - A(R - R 2 ) , (2.22) 
(HFB | A>(// c ) | HFB) = iV c (/i c ) = N. (2.23) 



Here E kin [ Pc ] = Ti(t Pc ) and E? nt [p c , v c ] = £f nt(normal) [p c ] + E c anomal [p c , z> c ] with 

K^JPc,i>c] = \(i>l^\Pc]i>c). (2-24) 



The above consideration has been carried out without any pre-assumptions 
concerning the density functional E int [R] E^ nt [ Pc , u c ] . Now, recalling the 
Hohenberg-Kohn theorem [27], we specify that the density functional can be 
chosen to be of a local form, i.e. dependent on the normal and anomalous local 
real densities p(r, r) and z/(r, r) defined in Appendix A by (A. 61) and (A. 64), 
respectively. Then -E r c lormal is a functional of the normal densities (isoscalar, 
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isovector, spin-orbital, etc.), and -Anomal is a functional of the normal and 
anomalous densities, the latter acquiring the simple form 

^nomaltPc^c] = £ / df U* c (f, r)^ (f, T; [ Pc ])u c (f, t) . (2.25) 
r=n,p J 



The pairing potential is defined by 

A(f, T) = l ^omajPc^c] (2 2g) 

2 5is c (r,T) 



That corresponds to the gap equation which takes on now a very simple mul- 
tiplicative form: 

A(f, r) = J*(f, r; \p c ])u c (f, r) . (2.27) 



The normal mean-field potential is given by 



U^r) = 5 ff- l f . (2.28) 



Due to the density dependence of it includes also a contribution arising 
from the variation of the pairing interaction energy: 

TJ . (jf T \ _ ^anomal [Pc,^c] ,„ ^ 



We emphasize that having found the gap function A(r), and the mean-field 
potential U(f), the Gor'kov equations [44] may be solved, for spherical nu- 
clei, exactly by using the coordinate-space technique (see Refs. [25,20,26] and 
Appendix C). The generalized Green's function obtained this way can be inte- 
grated over energy up to e c to yield both the normal and anomalous densities 
p c and u c which are used then to compute the energy of the system. This 
corresponds to a full HFB treatment of the nuclear ground state properties. 
It remains, of course, an art to find a "good" local functional, in particular its 
anomalous part. 



3 The normal part of the functional 



In this and the next section we describe in some detail the energy-density 
functional which we use in the present calculations. We shall omit in the 



12 



following the cutoff index c for simplicity. The functional is local in the sense 
that it contains only normal and anomalous densities, not the density matrix. 
The interaction part of the energy density is 

£int = c main + £Coul + £sl + £ss + c-anomal • 

(3.1) 



The terms e ma _ in and £coui have been described in Refs. [21,45], where also 
the connection of the parameters involved with the characteristics of nuclear 
matter has been given. Two different versions of the spin-orbit part e s i have 
been used: the one of Ref. [21], and also, after [20], the variantQ 

e sl {r)= l -C rl £ L ik [Vp { • ^(r, s)\p x ^(r, 

i,fc=p,n I ss' 
ss' 

+ 9i EE^V. sKwHr, s')Y E(^ f (r, s)a aP ^(r, s')) k 1 (3.2) 

a,/3 ss' ss' ) 

where ijy and ip are field operators, the brackets denote the ground state 
expectation value, and a superscript k indicates that a projection on particles 
of type k is performed. 

As given by (3.2), e s i corresponds to the two-body spin-orbit interaction 

F$ = C r 2 { K + K'f V T 2 )[V x 5{fi - f 2 ) x [p x -^)]-(o ? i + <? 2 ) , (3.3) 



and to the first-order velocity harmonic of the spin-dependent part of the 
Landau interaction amplitude 

= CVo(#i + g'iTi-T 2 )5(fi - f 2 )(ov<? 2 )(p r p 2 ) , (3.4) 



which has to be symmetrized in such a way that the momentum operators pi |2 
never act on the delta-function 5{f\ — r 2 ). 

In these expressions, the multiplier r 2 , has been introduced to make the param- 
eters k and gi dimensionless. It is given by r 2 . = (3/87rp ) 2 / 3 where p is the 
equilibrium density of one kind of particles in symmetric nuclear matter. The 
factor C is the inverse density of states at the Fermi energy, e F = k 2 . F /2m*, 
in saturated nuclear matter. It is given by C = 2e p/3p = 7r 2 / k 0F m* . 

7 As a rule, the convention H = 1 is used throughout the paper but fi will appear 
in some expressions. Hopefully this should not lead to any confusion. 
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Defining spin-orbit densities by 



then in spherically symmetric nuclei, the spin-orbit energy density (3.2) can 
be simplified to [20] 

= C ^o i E p (lp> lk ^T + -tz&tfJk) • (3-5) 

In deformed nuclei, the expression for e s i cannot be written in such a sim- 
ple form. In the axially symmetric case the derivation can be performed in 
cylindrical coordinates; the resulting formulae are given in [23]. 

We just mention some features of those parts of e- m t which are not given here in 
detail: The main contribution £ ma j n consists of a volume part and of a surface 
term both containing simple fractional- linear functions of the normal isoscalar 
density p + = p n + p p . This leads to the effective density-dependent forces in 
the scalar-isoscalar and scalar-isovector channels (see e.g. Refs. [20,21,45] for 
more details). There is no momentum dependence in these parts, therefore 
we have the effective mass m* equal to the bare mass m. The surface term 
involves a finite, but small Yukawa range parameter R = 0.35 fm; in the limit 
R — > it could be written with the help of the Laplace operator acting on 
functions of the density. It vanishes for constant density. The Coulomb part 
ecoui is taken in the usual form [19] with an exchange contribution in the Slater 
approximation. 

Note that in the more complicated parts of the force like those of eqs. (3.3) 
and (3.4) we did not introduce any density dependence. In general, all of the 
constants k, g±, and others to follow below could be density dependent as well. 
Until now there are no experimental data which would make it necessary to 
introduce such a dependence for the laboratory nuclei within or not too far 
from the stability valley. However, going to the nucleon drip lines and beyond, 
a more complicated density-dependent spin-orbit force might be of relevance. 
This is indicated by the relativistic mean field models for very neutron-rich 
nuclei [46] and by the exact Monte Carlo methods for small pure neutron 
drops [47] . Both approaches predict a much smaller (by a factor of 2-3) spin- 
orbit potential than the usual Skyrme-type ansatz. The role played by the 
density dependence of the effective spin-orbit interaction is left to be studied 
in future work. 

The spin-spin term e ss has hitherto been omitted since its expectation value in 
spin-zero nuclei vanishes. We include it here for completeness; it might become 
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important for magnetic excitations or polarization. A possible simple form of 
it is 



Note that the strength oc can be considered as the zeroth term in a decom- 
position of the spin-spin interaction amplitude at the Fermi surface in Legen- 
dre polynomials of the scattering angle. The next (first) term, proportional to 
g\ k , has been included in the spin-orbit part of the functional; however, the 
normalization chosen in eq. (3.2) is different from that of Ref. [10]. 



As long as there is no static spin polarization ((a) = 0) in the ground states 
of even nuclei and our interest is only in electric and mass quantities we can 
continue to neglect e ss . One delicate point should be addressed in this con- 
nection. Usually, when deriving the Hartree-Fock equations, one starts with a 
Hamiltonian containing a well behaved and properly antisymmetrized inter- 
action. In the case of contact forces, antisymmetry prevents self-interaction 
of the particles, so there is no need for precautions against that. Going over 
to an effective force which, as in our case, is not antisymmetrized in the ph 
channel, the self-interaction is not automatically excluded. Dealing with a 
not-too-small number of particles, by fitting to the experimental data the self- 
interaction is compensated (or accounted for) in the average by the choice of 
the parameters. 



This does not hold for the spin-spin part of the force. For an odd nucleus, 
consisting of a core with J = plus one particle, there is a nonvanishing spin 
density, and from eq. (3.6) there is also a contribution to the total energy 
which to lowest order is just the self-interaction of the odd particle and should 
therefore be excluded. On the other hand, considering the contribution of this 
part of the functional to the single particle potential when determining the 
wave function of the core, will give the polarization of the core which is a 
physical quantity. - - Usually, on the level of Hartree-Fock type calculations 
the spin-spin interaction is left out; it may be considered in a second step by 
perturbation methods. 





(3.6) 



i,fc=n,p 
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The basic set of functional parameters in the notations of Refs. [20,21,48] used 
in the present calculations is the set DF3: 



a\ = -6.422, 

h\ + = 0.163, h v 2+ = 0.724, 

a% = -11.1, h% = 0.31, 

r 



i>p = 0.285, n pn = 0.135, 
= 1.147 fm, 



a v _ = 5.417, 

h\_ = 0, h\_ = 3.0, 

a s _ = -4.00, h s _ = 0, 

9T = -g? = -o.i2, 

R = 0.35 fm. 



(3.7) 



This set set has been specifically deduced in [48] to reproduce not only the 
already known ground-state properties of magic nuclei but also the very recent 
experimental data [49] on the single-particle energies near the "magic cross" 
at 132 Sn. The set (3.7) corresponds to the following nuclear matter character- 
istics at saturation: compression modulus K = 200 MeV, chemical potential 
(binding energy per nucleon) /i = —16.05 MeV, asymmetry energy parame- 
ter p = 28.7 MeV, saturation density 2p = 0.1582 fm" 3 (A; 0F = 1.328 fm" 1 , 
e 0F = 36.57 MeV, C = Tr 2 /k 0F m = 308.2 MeV-fm 3 ). As in Ref. [20], the func- 
tional contains zero-range isoscalar spin-orbit interaction oc k and velocity- 
dependent spin-isospin interaction oc g[ (i.e. the first Landau-Migdal harmonic 
in the ar channel). This is in contrast to [21] where a finite range spin-orbit 
force had been used. 

With the above parameter set of the normal part of the density functional the 
ground states of magic nuclei are described fairly well. For example, the calcu- 
lated rms charge radii {r 2 ch ) 1/2 of 40 Ca, 48 Ca and 208 Pb are 3.480 fm, 3.478 fm, 
and 5.500 fm, respectively, which agree nicely with experimental values de- 
duced very recently [50] from the combined analysis of optical, muonic and 
elastic electron scattering data: {r 2 cYl ) l/2 = 3.4767(8) fm for 40 Ca, 3.4736(8) fm 
for 48 Ca, and 5.5013(7) fm for 208 Pb. The predictions for some other double 
magic nuclei are: (r 2 ch ) l l 2 = 3.721 fm ( 56 Ni), 3.951 fm ( 78 Ni), 4.453 fm ( 100 Sn), 
and 4.705 fm ( 132 Sn). 



4 The pairing part of the functional 

The pairing energy density £ an omai in eq. (3.1) is chosen in the form prescribed 
by (2.25): 

^ a „o m ai(r1 = £ T^(f; [p})\v T {f)\ 2 . (4.1) 

r=n,p 
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As discussed in Refs. [4-6] and mentioned in the introduction, effective three- 
or more-body forces are indispensable if one wants to reproduce isotope shifts 
in charge radii which means that £ a nomai must depend on the normal den- 
sity [22,24], and therefore we introduce 

<F e,nn = <F e,p P = C7o ^ a .) j ( 42 ) 



The dimensionless strength is supposed to be, within the EDF approach, 
a local functional of the isoscalar dimensionless density x = (p n + p p )/2p , 
with p n ( p ) the neutron (proton) density; C and p are defined in the previous 
section. We shall use the parametrization of suggested in Ref. [24]: 

f(x) = f(x(r}) = /| x + h^(f) + f%rl (Vx(f)) 2 . (4.3) 



The parameter /| x is negative and simulates an attraction in the pp channel in 
the far nuclear exterior, and / v are taken to be positive [24] . The exponent 
q in the second term is introduced to have a more flexible parametrization. 
A repulsive short-range part of the G-matrix effective interaction with a x 2 ^ 3 
dependence has been discussed, e.g., by Bethe many years ago [51] while de- 
veloping a Thomas-Fermi theory for large finite nuclei. The choice q = | 
seems thus to be reasonable, and indeed our calculations showed that some 
improvements in reproducing the isotopic shifts may be achieved with this 
choice compared to the linear case q = 1. We shall use q — | in the present 
paper. The three-body force of Refs. [5,6] corresponds to a linear dependence 
of J* on p (q = 1, /| = 0). As shown in Ref. [24], the self-consistent EDF 
(HF+BCS) calculations with the density-gradient term oc /y in pairing force 
provide desirable size of isotopic shifts and right order of odd-even staggering 
observed in lead isotopes, the coupling of the proton mean field with neutron 
pairing being the major effect in this case. The variation of the anomalous 
energy (4.1) incorporating the above force with respect to normal densities 
gives the following contribution to the central mean-field potential: 



- — ^qh^x q - 1 - 2r 2 4 (Ax + Vx ■ v)] (k| 2 + \u v \ 2 ) . (4.4) 



Upair 4p 



As will be shown in the next section, different choices of the parameters of 
the particle-particle force (4.3) are possible to reasonably describe the neu- 
tron separation energies and isotopic shifts in charge radii. In particular, the 
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following sets are deduced for the lead isotopes: 



J ex 


-0.56, 


7 C 

h> = 


0. 


r-f 
/v 


= 


(a) 


f ? = 

J ex 


-1.20, 


ht = 


0.56, 


/v 


= 2.4 


(b) 


f ? = 

J ex 


-1.60, 


tf = 


1.10, 


/I 


= 2.0 


(c) 


f ? = 
J ex 


-1.79, 


tf = 


1.36, 


/J 


= 2.0 


(d) 


f ? = 
J ex 


-2.00, 


ht = 


1.62, 


/I 


= 2.0 


(e) 


= 

J ex 


-2.40, 


ht = 


2.16, 


/I 


= 2.0 


0, 



(4.5) 



Let us, in the rest of this section, study the behavior of A as a function of 
density p = 2p x (or the Fermi momentum kp = (37r 2 p/2) 1//3 = k^x 1 ^ 3 ) hi 
symmetric (p n = p p ) uniform infinite nuclear matter with pairing interaction 
of eqs. (4.2), (4.3), and also discuss how the pairing affects the equation of 
state (EOS) and the position of the equilibrium point. The density-gradient 
term oc /y vanishes in this case, thus only the terms with parameters /| x and 
ht are left. For uniform matter the gap equation (2.27) reduces to 

a / \ ^t, ^ f d/c A(x) , „ „. 

A x = -J*(x) / K J , 4.6 

Jk< kc (27r) 3 2v /( efc - eF (x)) 2 + A 2 (x) 



whereQ k c = J 2m(eF + e c )/h and = h 2 k 2 /2m. Canceling A(x) on both 
sides of this equation, to find a nontrivial solution A^0, one gets 

-x^Hx) )\t ^ =-1, (4.7) 



where 8{x) = A(x)/ep(x) and t c (x) 
measured from the Fermi energy e-p : 



= 1 + e c /ep(x) with e c the energy cutoff 
€ fx 2j/3 . The solution of the last equation 



in the weak pairing approximation, for f^(x) < 0, is given by (see, e.g., [25,52]): 



A(x) = 8cof x 



2/3 



S(X) 



\ s(x) + 1 



exp s(x) — 2 + 



ft(x)x 1 / 3 / 



(4- 



8 Note that, in our approach, the upper limit in the truncated momentum or energy 
space depends on density since e c is chosen to be measured from the Fermi level 
position determined by the chemical potential \x (see Appendix C). At any given 
density, p is introduced as a Lagrange multiplier and therefore, in varying the cutoff 
functional, the phase space is kept fixed. 
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with s(x) = y 1 + e c /e-p(x) = k c /kp. In the far nuclear exterior when x — > 0, 
JF^ should be determined by the free NN scattering. However, approaching 
the nuclear surface, finite range and nonlocal in-medium effects may already 
considerably modify the force even at low density, and therefore the best 
parametrization of an effective contact force at 2; < 1 might not necessar- 
ily agree with the vacuum value which could be extracted from the free NN 
interaction, as discussed by Migdal many years ago [10]. 

Shown in Fig. 1 are the values of dimensionless strength ft (upper panel) 
and the results for A (lower panel) in infinite matter with the parameter sets 
(a)-(f) of eq. (4.5) deduced for the lead chain. The calculations are performed 
with cutoff e c = 40 MeV and Fermi energy cof = 36.57 MeV of saturated 
nuclear matter from the parametrization (3.7) of the normal part of the density 
functional. 

One should notice that eF entering the integrand of the gap equation (4.6) 
can be expressed directly through density by ep = fi hp /2m with kp = 
(3tt 2 p/2) 1/3 = k 0F x 1/3 only if the pairing is weak indeed so that the dependence 
of the Fermi energy (and the chemical potential /i) on A can be disregarded. 
Otherwise one should introduce the particle number condition 

2 f dk . . 
p Jk<k c {2-Ky i 



where 



and solve the system of the two equations (4.6) and (4.9) with respect to A 
and eF- The results shown in Fig. 1 by full lines correspond to such a direct 
solution while those shown by dashed lines to the weak pairing approximation, 
eq. (4.8). 

As seen in Fig. 1, the approximation (4.8) works well in the entire range of 
kp for the set (a), but for the other sets this is true only at kp greater than 
«1.2 fm _1 and also, for the sets (b), (c) and (d), at kp less than 0.42, 0.14 
and 0.042 fm -1 , respectively (in these regions the ratio A/ep does not exceed 
0.1). For the sets (e) and (f), at lower densities, eq. (4.8) can not be used any 
more to estimate the gap even by order of magnitude since, as seen from the 
behavior of the corresponding dashed lines, it becomes divergent. 

All parameter sets (4.5) except (a) reproduce the neutron separation energies 
and the isotope shifts of charge radii of lead isotopes reasonably well (see next 
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section). For the sake of comparison also shown in Fig. 1 are the values of the 
l So pairing gap in nuclear matter obtained for the CD-Bonn potential with- 
out medium effects (using free single-particle spectrum e& = k 2 /2m) [53] and 
for the Gogny Dl force in the HFB framework [54]. The agreement between 
the two latter calculations is relatively good while both deviate noticeably 
from our predictions. The curve for contact density-independent pairing force, 
set (a), stands by itself with a positive derivative dA(x)/ da; everywhere; no 
acceptable description of (r 2 h ) could be obtained in this case (see Fig. 5). 




Fig. 1. Dimensionless pairing strength (upper panel) and pairing gap (lower panel) 
in infinite nuclear matter as functions of the Fermi momentum. In the upper panel, 
the horizontal arrow shows the critical value of /| x (/| r = —1.912) at which an / = 
spin singlet nucleon-nucleon bound state appears at zero energy while the vertical 
one marks the Fermi momentum at saturation point for the functional DF3 without 
pairing, /cof = 1.328 fm _1 . In the lower panel, full (dashed) curves (a)-(f) are 
obtained by solving eqs. (4.6) and (4.9) (by using the weak pairing approximation, 
eq. (4.8)) with contact pairing force (4.2), and an energy cutoff e c = 40 MeV, 
and correspond, respectively, to the parameter sets (a)-(f) of eq. (4.5). The solid 
circles and stars are the solutions of the nonlocal gap equation with the CD-Bonn 
potential [53] and with the finite-range Gogny Dl force [54], respectively. The dotted 
line shows the values of the pairing gap calculated using eq. (4.22) with the free NN 
scattering phase shifts (see text). 



20 



Interestingly, for the sets (b)-(e) which reproduce satisfactorily both the neu- 
tron separation energies S n and mean squared charge radii (r^ h ) there exists 
a pivoting point at k F fa 1.15 fm _1 (at ^0.65 of the equilibrium density) with 
the same value of A piv « 3.3 MeV (/| iv « -0.78). 

One can see in Fig. 1 that for parameter sets with bigger absolute values of 
and the region where A(x) varies strongly moves towards lower densities, 
the slope becomes steeper and around x ~ 1 the pairing gap tends to be 
very small indeed. One may think that in finite nuclei the effect of the strong 
dependence of A on x, especially at small densities, would not influence the 
nuclear properties noticeably. But even if A would be concentrated only in 
the surface region and outside of the nucleus, the anomalous density, due to 
quantum effects, would not vanish inside the nucleus where x ~ 1. In this 
case the effect of density dependence should persist in the nuclear volume and 
should be more pronounced at larger values of ht. This point was confirmed 
by our calculations (cf. Ref. [24] and below). 

Consider in more detail the behavior of A at very low densities when kp — > 
assuming the weak pairing regime, i.e. /| x > /| r where /| r is a critical strength 
constant given by 

/£ = -2^, (4.H) 



with k 0c = ^2me c /h. This constant is determined by the condition 




(4.12) 



obtained from the gap equation (4.6) at A = e F = 0. For the chosen energy cut- 
off and parametrization (3.7) of our density functional we have /| r = —1.912. 
To leading order, on the other hand, at kp — > from eq. (4.8) we obtain 

A = C£FeXp (2^)' a<0 ' (413) 



where c = 8e 2 ~ 1.083 and where we have introduced the quantity 

7r / v / 2me c 

a = 

2/c F V ^OF 





(4.14) 



It can be shown that this quantity is nothing but the singlet nucleon-nucleon 
scattering length. To be more specific, consider first the two-neutron problem 
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in the vicinity of the critical point, /| x ~ /| r , when the attraction is strong 
enough to produce a bound pair state (the scattering problem at any ft will 
be considered below). The bound-state wave function for the relative motion 
of two neutrons in case of contact interaction C f^ x 5(r) is determined by (cf. 
Ref. [55]): 

^(f) = -C fLG (r, e b )V(0) , (4.15) 



where G is the free Green's function in the truncated space, 

dk e y 



Go(r,e h )= -f-^ , (4.16) 

Je™<er (2vr) 3 e£ n -e b 



with e b the binding energy and e£ n = h 2 k 2 /m (the reduced mass is ra/2). In 
the rest system of a nucleus, the center-of-mass energy e nn in the scattering 
problem would correspond to the energy e nn /2 of each of the two nucleons 
in the s-wave pairing problem. This implies that the cutoff in the /c-space in 
eq. (4.16), /c° n = ^/me™//j, should be the same as in the gap equation. Thus, 
in the energy space, one has e° n = 2e c and from (4.15) at r = one obtains 
the equation to determine e b : 



1 = -^/ de ^b- < 4 - 17 > 



V^F cx { e-e b /2' 



which reduces to 



^=( i -i/?«*»i/5)"- < 4i8) 



where /| r nn is a critical value of /| x at which eq. (4.15) has a bound state 
solution e b = 0. Comparing eq. (4.17) at e b = with eq. (4.12) one finds 
that /| r nn coincides with the value defined above by eq. (4.11). Now, in the 
vicinity of /| r we can set arctan ^/2e c /|e b | tt/2 and use e b = —h 2 /ma 2 m in 
eq. (4.18). Then it is easy to see that the solution for the scattering length 
a nn is exactly the same as given by (4.14). The expression (4.13) for A agrees 
with the results of Ref. [40] based on a general analysis of the gap equation at 
low densities when fc F |a| I. But we should stress that (4.13) is valid only 
in the weak-coupling regime corresponding to negative a. In the opposite case 
the gap in the dilute limit has to be found in a different way. 

At /| x > /| r the scattering length is negative, and from eq. (4.13) it follows 
that at low densities the pairing gap is exponentially small and eventually 
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A(kp — > 0) = 0. Such a weak pairing regime with Cooper pairs forming 
in a spin singlet I = state exists up to the critical point at which the 
attraction becomes strong enough to change the sign of the scattering length. 
Then the strong pairing regime sets in, eqs. (4.8) and (4.13) are not valid any 
more, and A should be determined directly from the combined solution of 
the gap equation (4.6) and the particle number condition (4.9). In the dilute 
systems, eF plays the role of the chemical potential [i. The latter is defined by 
\i = €f(A;f) + U(kp) with U(k-p) the HF mean field at the Fermi surface which 
is negligible for the fermion gas. At the critical point \x becomes negative and 
a bound state of a single pair of nucleons with the binding energy eb = 2/i 
becomes possible [56,57]. This can be easily seen from the gap equation (4.6) 
written in the form 

(k 2 \ i r dk' 

- - 2„) * = -sg„( £t - -^J** , (4.19) 



where we have introduced the functions (f>k = A/y (e& — /i) 2 + A 2 and replaced 
€f by /i. In the strong coupling regime, /i < 0, and in the dilute limit, \<f>k\ <C 1, 
this equation reduces to the Schrodinger equation for a single bound pair 
where 2/x plays the role of the eigenvalue. It is equivalent to the coordinate- 
space equation (4.15). Thus 2/x = eb where the binding energy eb is determined 
from (4.18). The pairing gap A at low densities in this regime can be found 
from the particle number condition (4.9) which in the leading order now reads 



In the vicinity of the critical point from this equation we find 

A 2 - — xe 2 1 

37r 0F k 0F a ' 



which gives 

\^y\ «>o. (, 21 , 



a = ^/2ttp\ 1 / 2 



From this consideration it follows that, in the dilute case, the energy needed to 
break a condensed pair goes smoothly from 2A to 2/j = e b as a function of the 
coupling strength when the regime changes from weak to strong pairing. But 
as seen from (4.13) and (4.21), the behavior of A at low densities is such that 
the derivative dA/ dp at p — > as a function of /| x exhibits a discontinuity 
from to oo. This is illustrated in Fig. 2 where we have plotted the gap 
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A(p) at very low densities for the sets (c)-(e) of eq. (4.5) embracing both 
regimes. We note also that the analytical expressions (4.13) and (4.21) give a 
purely imaginary gap at the critical point when the scattering length changes 
sign {k F a — > y/2mjl/h — > i if /i becomes negative). In this connection it is 
instructive to write the weak coupling expression (4.8) for A in the following 
form [58]: 



A (/cp) = ce F exp 



TT 



cot 5 (h?) 



(4.22) 



where c = 8e 2 and where we have introduced the Fermi level phase shift 
6 (hp) defined by 

/ +xn \ 4k ° F ( 1 ,k c (k F )\ k F fk c (k F )-k F \ fAO „. 

k F cot S(k F ) = + — In , x , , , 4-23 

n \ft{k F ) 2k 0F J n \k c (k F ) + k F J 



with k c (k F ) = y /c 2 c + k F . Eq. (4.23) corresponds to an exact solution of the nn 
scattering problem at the relative momentum k = k F with the states truncated 
by a momentum cutoff k c = k c {k F ) for contact interaction Cof^(k F )5(r) (see, 
e.g., Ref. [59]). In the dilute limit, from eq. (4.22) one notices again that the 
pairing gap becomes pure imaginary at the critical point since cot 5 i when 
|a nn | — > oo. At very low densities, with the parametrization (4.3) of the pairing 
force and with the chosen density-dependent cutoff, eq. (4.23) reduces to 




, (4.24) 



P (fm 3 ) 

Fig. 2. Pairing gap A as a function of p at very low densities. Curves (c)-(e) corre- 
spond to the parameter sets (c)-(e) of eq. (4.5), respectively. 



24 



where a nn is the scattering length defined by eq. (4.14) and r nn is the effective 
range, r nn = 4/irk 0c . The first two terms in this equation would describe low- 
energy behavior of the nn s-wave phase shift through an expansion of k cot 5 
in powers of the relative momentum k = hp if the interaction were density- 
independent - in our case, if the coupling strength and momentum cutoff were 
fixed by ft = /| x and k c = k Qc , respectively. It follows that with a density- 
dependent effective force, such an expansion contains additional terms which 
for the parametrization used here are of the same order as the effective range 
terrn^]. This simply demonstrates that, for reproducing the pairing gap, the 
effective interaction even at very low densities need not necessarily coincide 
with the bare NN interaction. 

When the Fermi momentum k F approaches from below the upper critical point 
at which the pairing gap closes, k F = k p(—f^ x /h^) 1 ^ 3g defined by f^(kp) = 0, 
we get from eq. (4.22) 

A(M^ce F exp - f ™ . (4.25) 

3qfS x {k F - k F r ) . 



Thus, at higher densities, the pairing gap becomes exponentially small when 
kp approaches k F , in agreement with general analysis of the gap solutions [40]. 
In weak coupling, as we have already shown, A(fcp) is also exponentially small 
at low densities, in the vicinity of kp = 0. It is noteworthy that near both 
critical points the pairing potential A can be found from the behavior of the 
phase shift by using eq. (4.22) with a smooth density- dependent prefactor 
c(/cf) (the details will be given elsewhere [61]). In the present paper, as an 
illustration, we show in Fig. 1 by the dotted line the values of A(fcp) obtained 
from eq. (4.22), with a prefactor c = 8e~ 2 , by using "experimental" nn phase 
shifts, without electromagnetic effectsp ]. It is seen that A obtained this way at 
low densities closely follows the solution of the gap equation with the CD-Bonn 
potential [53]. The nn phase shift passes zero at the relative momentum k w 
1.71 fm -1 , and the gap should vanish at the corresponding Fermi momentum. 
Unfortunately, the solutions for A are given in Ref. [53] only in the region up 
to kp = 1.4 fm -1 which is rather far from the upper critical point. 



9 Our effective pairing interaction with the choice q = 1/3 would lead in the dilute 
limit to the expression for A of the form of eq. (4.13) but with a different prefactor c 
depending on the value of h^. If, furthermore, we define the latter parameter by = 
(l + 21n2)(/! x ) 2 /6 we get in the leading order A(k F ) = (2/e) 7/3 e F exp(vr/2A; F a), i.e. 
the result obtained in Ref. [60] for a non-ideal Fermi gas by studying the singularities 
of the interaction amplitude (vertex function T) and taking into account the terms 
up to the second order in k F a. 

10 We thank Rupert Machleidt for providing us with these nn phase shifts. 
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For symmetric nuclear matter, with the functional DF3 used in the present 
paper, the energy per particle is 

E . . 2 r dk Ti 2 k 2 . . 1 v i'\ / \ SA 2 (x) . . 

- (x) = / n k {x) + -e 0F a v + fl(x)x + V , 4.26 

A p x Jk<k c (2ny 2m 3 2j^(x)xeoF 



where /+(#) = (1 — /(l + /i2+ x ) with the parameters given by (3.7); here 
the "particle-hole" term oc vanishes in the dilute limit linearly in density. 
The chemical potential is 

p(x) = e F (x) + -e 0F al[fr(x)x 2 + 2fl(x)x] + ^ { J J > , (4.27) 



where the prime denotes the derivative with respect to the dimensionless den- 
sity x. The Fermi energy e F (x) and the pairing gap A(x) entering these equa- 
tions are determined from (4.6) and (4.9). The last two terms in (4.27), even 
in strong pairing regime, vanish in the dilute limit at least as x q if < q < 1 
or linearly in x if q > 1 (in our case, the exponent in the density-dependent 
pairing force is q — |). Thus, we see again that, in strong coupling, in the 
leading order \i = e F = £h/2 < 0. 

The calculated energy per nucleon as a function of the isoscalar density p is 
shown in the upper panel in Fig. 3 together with the results of the nuclear 
matter calculations [62,63] for the UV14 plus TNI model. It is seen that DF3 
gives qualitatively reasonable description of the nuclear matter EOS and that 
pairing could contribute noticeably to the binding energy especially at lower 
densities. In the lower panel in Fig. 3 we have plotted the pairing energy per 
nucleon, (E/A) pSuir , obtained by subtracting from eq. (4.26) the corresponding 
value of E/A at A = 0. This pairing energy is a sum of the positive contri- 
bution coming from the kinetic energy, i.e. from the first term in (4.26), and 
the negative anomalous energy - the last term in (4.26) (an expression for 
(E/A) pair in the case of weak pairing is derived in Appendix B). The pairing 
contribution is small for the density-independent force with = 0.56, set (a) 
of eq. (4.5), and increases, as expected, for the sets (b)-(f) with a shift to 
lower densities as /| x becomes gradually more attractive. For the sets (e) and 
(f) the attraction is strong, /| x < /| r . In these cases a nonvanishing binding 
energy in the dilute limit is solely due to Bose-Einstein condensation of the 
bound pairs, the spin-zero bosons, when all the three quantities, p, E/A and 
(E/A) piliT , reach the same value e^/2 (e& = —0.0646 and —1.616 MeV for the 
set (e) and (f), respectively). This is illustrated in Fig. 4 where we have plotted 
E/A and (E/A) pair as functions of p at very low densities. Analytically, for 
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Fig. 3. Energy per nucleoli E/A (top) and pairing contribution to E/A (bottom) 
in symmetric nuclear matter. Curves (a)-(e) are calculated using eq. (4.26) and 
correspond to the strength parameters (a)-(e) of eq. (4.5), respectively. Open circles 
and stars are the calculations of Ref. [62] and Ref. [63], respectively, for the UV14 
plus TNI model. The vertical arrows mark the saturation density for the functional 
DF3 without pairing, 2po = 0.1582 fm -3 . 

the kinetic energy term in (4.26) in the leading order we find 
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where we have used the definition (4.11) and the relation 
Combining this with the last term in (4.26), one gets 
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By using eq. (4.21), this expression reduces exactly to —% 2 /2ma 2 = It 
follows that the model just described is in fact parameter-free: in the strong 
coupling regime near the critical point, the ground state properties of nuclear 
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Fig. 4. Energy per nucleon E/A (top) and pairing contribution to E/A (bottom) for 
symmetric nuclear matter at low densities. The notations are the same as in Fig. 3. 

matter in the dilute limit are completely determined by the scattering length. 

We have considered subsaturated nuclear matter with 1 Sq pairing within our 
local EDF framework and demonstrated some results, including extrapolation 
to the low density limit, with a few possible parameter sets of the pairing 
force deduced from experimental data for lead isotopes. At low densities, in 
the T = case of symmetric N = Z matter, the 3 Si — 3 D\ pairing corre- 
lations leading to a Bose deuteron gas formation might be more important 
since the n-p force is more attractive than that in the p-p or n-n pairing 
channels (see Ref. [64] and references therein). Thus, our approach, with 1 S'o 
pairing only, would be more relevant for an asymmetric N ^ Z case and for 
pure neutron systems. From this point of view the best choice of the effec- 
tive contact pairing force for the EDF calculations seems to be the set (d) 
of eq. (4.5) since it gives the singlet scattering length a nn w —17.2 fm which 
corresponds to a virtual state at ~ 140 keV known experimentally^]. As seen 

11 The contact interaction in a truncated space of states, with two parameters /| x 
and fcoc 

can be calibrated to produce not only the empirical scattering length ct nn 
but also a realistic effective range r nn ~ 2.8 fm which is directly related to the 
momentum space cutoff by r nn = 4/irko c , see Ref. [59]. It would yield a rather low 
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in Fig. 1, with this choice the behavior of A at low densities agrees well with 
the calculations based on realistic NN forces. At higher densities, however, 
our predictions for A with the set (d) go much higher reaching a maximum 
of 4.84 MeV at k F ps 0.92 fm _1 while the calculations of Ref. [53] give a 
maximum of about 3 MeV at kp ps 0.82 fm _1 . With a bare NN interaction, as- 
suming charge independence and m n = m p in the free single-particle energies, 
the pairing gap would be, at given kp, exactly the same both in symmetric 
nuclear matter and in neutron matter. As shown in Refs. [65,66], if one in- 
cludes medium effects in the effective pairing interaction, most importantly 
the polarization RPA diagrams in the cross channel, the pairing gap in neutron 
matter would be substantially reduced to values of the order of 1 MeV at the 
most. Whether such a mechanism works in the same direction for symmetric 
nuclear matter is still an open question. If the gap were smaller than that 
obtained with the Gogny Dl force, which is also shown in Fig. 1, it would 
be difficult to explain the observed nuclear pairing properties. The effective 
contact density-dependent force (4.3) with our preferable parameter set (d) of 
eq. (4.5) yields larger pairing energy in nuclear matter than the Gogny force, 
but in finite nuclei this is compensated by the repulsive gradient term. The 
phenomenological pairing force used here contains dependence on the isoscalar 
density only since we have analyzed the existing data on separation energies 
and charge radii for finite nuclei with a relatively small asymmetry charac- 
terized by (N — Z)/A < 0.25. An extrapolation to neutron matter with such 
a simple force would give a larger pairing gap than for nuclear matter. This 
suggests that some additional dependence on the isovector density p n — p p 
might be present in the effective pairing force. This possibility is planned to 
be tested in our future work, with a more careful analysis of experimental data 
though the relevant data base is not rich enough. 

Now consider how the density changes when the pairing gap appears in nuclear 
matter. To leading order, one finds the following expression for the energy per 
nucleon near the saturation point: 



where K is the compression modulus at saturation density, /3(x)I 2 is the 
asymmetry energy, with / = (p n — p p )/2p x = (N — Z)/A. This expression 
is valid at \x — x \ <C 1, |/| <C 1 and A C e F . Higher order effects connected 
with /-dependence of K, j3, A and ep are neglected. The derivation of the 

value for the energy cutoff e c ps 9 MeV. However, as we have already mentioned, 
even at low densities the force can be considerably modified by nonlocality and 
polarization effects so that the energy dependence of the Fermi-level s-wave nn 
phase shift, eq. (4.24), at low relative momenta k governed by the effective range 
might be different from the free two-nucleon case. 
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last term in (4.30) is discussed in detail in Appendix B. Due to this pairing 
term, the position of the equilibrium point may be shifted to lower or higher 
densities depending on the behavior of A (a;) near x = 1. If A does not depend 
on x in the vicinity of x — 1 then the equilibrium density decreases due to 
presence of e-p{x) oc x 2 ^ 3 in the denominator of the pairing term (the system 
gains more binding energy). This means an expansion of the system. The 
effect is enhanced if A becomes larger during such an expansion, i.e. when the 
derivative dA(x)/ dx at x = 1 is negative. This point may be illustrated by 
the following simple consideration. At equilibrium the pressure P = which 
means 

£-(E/A)=0. (4.31) 



For saturated symmetric nuclear matter without pairing the dimensionless 
density is, by definition, xq = 1. When 1^0 and A^O, from eq. (4.30) with 
condition (4.31) the new equilibrium density can be found (in units of 2p ), 
which is the solution of the equation 
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(4.32) 



From this equation one can see that the density should be sensitive indeed to 
the derivatives of A, and a negative slope in A should cause a decrease of the 
density. The obtained relations can be used to estimate the influence of pairing 
on the charge radii for heavy nuclei as was done in Ref. [22]. It was shown 
that pairing interaction with strong p-dependence at x ~ 1 might significantly 
change the equilibrium density. For the parametrization used here the size 
of this effect is controlled by the parameter h$. The shift of the saturation 
point is relatively small: \8p\/2p < 0.8 % for all the six sets of eq. (4.5) 
(8p/2p ps +0.4 % and ps —0.4 % for the set (a) and (d), respectively). In 
finite nuclei, the surface term oc / v is equally important to produce a kink 
in the radius evolution along isotope chain at magic neutron number and, 
especially, to explain the observed odd-even staggering in (r^ h ) [24]. Although, 
as seen in Fig. 1, the calculations with bare NN interaction or with Gogny 
force give a negative sign for dA(x)/ dx near x ~ 1 (kp ~ 1.33 fm _1 ), but 
the slope might be not steep enough. This is the probable reason why the 
HFB calculations with the Gogny force, which give a good description of 
the global pairing properties of nuclei, could not reproduce the kink in lead 
isotopes [67]. The density dependence of the pairing force leads to the direct 
coupling between the neutron anomalous density and the proton mean field as 
given by (4.4). The suppression of |z/ n | 2 in the odd neutron subsystem because 
of the blocking effect influences the potential t/p air of the proton subsystem 
through the volume, oc ht, and surface, oc /y, couplings, and this moves the 
behavior of the proton radii towards the desired regime [24]. 
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5 Numerical results for some isotope chains 

The calculations for finite spherical nuclei are performed using the coordinate- 
space technique which is described in detail in Appendix C. In the con- 
struction of the Gor'kov Green's functions, the four linear-independent so- 
lutions (ui,Vi) ,i = 1-4, were found by the Numerov method with a radial 
step of 0.1 fm and physical boundary conditions imposed for each Ij channel 
at the origin and at r = 25 fm. The contour in the complex energy plane 
with Ime = ±8 MeV along the horizontal sections and with energy cutoff 
e c = 40 MeV measured from the Fermi level was used (see Figs. 21, 22 in 
Appendix C). The integration along the contour was performed by the Simp- 
son method with automatic step selection. The convergence of the iteration 
procedure is controlled by a few criteria: for two successive steps % and % + 1, 
the conditions \p i+ i(r) - Pi(r)\ < 10~ 6 fm -3 and |A m (r) - Aj(r)| < 50 keV 
should be achieved for all r, the chemical potential \i should finally satisfy 
the condition \N(fi) — N\ < 0.01 and the contribution Nij to the total parti- 
cle number N(p) from the states with higher angular momenta Ij should not 
exceed 0.01. Such criteria guarantee an accuracy not worse than 0.1% for all 
the calculated quantities of our interest. The mean square charge radii were 
computed from ground state charge densities obtained by folding the point nu- 
cleon distributions with nucleon charge form factors, including the relativistic 
electromagnetic spin-orbit correction, in the same way as in Ref. [68]. 

Let us discuss first the results obtained for the lead chain. As shown in the 
upper panel of Fig. 5, the experimental neutron separation energies are re- 
produced equally well for all the chosen parameter sets of the pairing force. 
The curves (a)-(f) correspond to the sets marked with the same letters in 
eq. (4.5). The curve (a) is obtained for the "constant" pairing force, without 
density dependence. The other sets (b)-(f) differ from that by non-zero values 
of and fy parameters. As already mentioned in Section 4, for each /| x value 
it is possible to find such values of ht and /y that the pairing energy for a 
given nucleus is practically the same. These parameters are kept fixed for all 
isotopes of the lead chain. 

In the lower panel of Fig. 5 the calculated isotope shifts of mean squared charge 
radii, S(r^ h ), for the lead chain with respect to 208 Pb as a reference nucleus are 
presented. These curves except (a) do not differ from each other significantly 
and all of them reproduce qualitatively the kink and the average size of odd- 
even staggering. The curve (a) corresponds to a simple pairing delta-force; in 
this case the behavior of S{r^ h ) is rather smooth and neither kink nor staggering 
are reproduced. The variant (a) of eq. (4.5) was included just to give an 
example that without any density dependence in the contact pairing effective 
interaction it is also possible to describe the neutron separation energies but 
not the evolution of charge radii. With the parametrization of eq. (4.3) used 
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Fig. 5. Upper panel: neutron separation energies S n for lead isotopes. Lower panel: 
differences of mean squared charge radii 5(r^ h ) with respect to 208 Pb as reference 
nucleus; 72% of the corresponding liquid drop values (using tq = 1.1 fm) are sub- 
tracted to enhance the visibility of the small differences. The calculations are per- 
formed with the functional DF3 and the coordinate-space technique. Curves (a)-(f) 
connect the points obtained with parameter sets (a)-(f) of the pairing force (4.5), 
respectively. Experimental data for S n , including those derived from systematic 
trends, are from [69,70]; data for 5(r^ h } are from [71-73]. 



here for the pairing force, the size of the staggering and the kink in charge radii 
of Pb isotopes can be reproduced only if the "gradient" parameter / v is ~ 2. 
Without gradient term the results for radii would be in between the curves (a) 
and (b) shown in the lower panel of Fig. 5 (an example is given in the lower 
panel of Fig. 6). On the whole, as seen in Fig. 5, the set (d) yields a somewhat 
better description than the other sets, particularly for the lighter Pb isotopes. 
For this reason, and due to the fact that the set (d) corresponds to the correct 
value of the singlet scattering length in the dilute limit (see previous Section), 
this set will be taken as our preferable choice, and the corresponding functional 
will be referred to as DF3(d) in the following. 

In Fig. 6 we compare the DF3(d) predictions for Pb isotopes with the HFB re- 
sults obtained for the Gogny D1S force [76] and for two state-of-the-art models 
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Fig. 6. The same as in Fig. 5 but only for the set (d) of the pairing force in com- 
parison with the HFB calculations with Skyrme forces SkP [74] and SLy4 [8,75], 
and the Gogny D1S [76] force. Also shown in the lower panel by the dotted 
line are the 6(r^) values obtained without gradient term in the pairing force 
(/| x = -1.79, = 1.66, /| = 0). 

based on Skyrme density-dependent forces SkP [74]and SLy4 [75]. The calcu- 
lations with Gogny force for the odd isotopes are done in the blocking approxi- 
mation, and thus we can show the results both for even and odd Pb isotopes^]. 
The Skyrme-HFB calculations^ were done, unfortunately, without blocking, 
and therefore we show the radii only for even isotopes; the S n values for odd 
isotopes were obtained by correcting the no-blocking HFB energy by adding an 
average pairing gap, this was considered to be a very good approximation. The 
details of the HFB+SkP and HFB+SLy4 calculations can be found in [77] (see 
also the references therein), and here we only briefly discuss how the pairing 
correlations have been implemented in these models. The HFB+SkP calcula- 
tions have been done with the same SkP force in the ph and pp channels, and 



We thank Jean-Frangois Berger, Jacques Decharge and Sophie Peru for providing 
us with these results. 

13 We thank Jacek Dobaczewski for providing us with a numerical hie of these 
Skyrme-HFB calculations. 
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therefore in these calculations the pairing force has density dependence rr 1 / 6 . 
The energy cutoff has been chosen by including the quasi-particle states up 
to the energy equal to the depth of the effective single-particle potential [74]. 
This implies that the cutoff depends not only on the particle numbers N and 
Z but also on the quantum numbers j and /; for /=0 it varies in the range 
of about 40-50 MeV. The HFB+SLy4 calculations have been done with the 
density-dependent zero-range pairing force in the pp channel, and the SLy4 
force in the ph channel. The form of this pairing force is identical to the 
velocity-independent piece of the Skyrme interaction (see Ref. [78]) with pa- 
rameters V = t = -2488.9f3 MeVfm 3 and V 3 = f9990 MeVfm 3+1 / 6 , and it 
also has the x 1 ^ 6 density dependence. The prescription for the energy cutoff 
has been identical as in the HFB+SkP calculations. We remark that such a 
prescription is much different from that which we use here; with such a recipe 
it is not easy to construct the corresponding effective contact pairing force 
with a fixed phase-space cutoff, hence no simple extrapolation to the pairing 
properties of uniform nuclear matter. 

One can see in the upper panel of Fig. 6 that the HFB+D1S, HFB+SkP, 
SkLy4 and our DF3(d) calculations describe the neutron separation energies 
with more or less the same quality, some deviations are observed in the region 
above 207 Pb. For the lighter Pb isotopes, the Gogny force slightly overestimates 
the odd-even effect in S n while with DF3(d) it is slightly underestimated. 

In the lower panel of Fig. 6 it is seen that the Skyrme-HFB calculations with 
SkP and SLy4 forces give too small kink in radii. Calculations with the Gogny 
D1S force for S(r 2 h ) could not reproduce the kink either. The latter model 
yields sizeable staggering in radii, but the effect is too small, at least by a 
factor of two smaller in amplitude than experimentally observed; moreover, 
its sign in isotopes below A=19A is reversed with respect to the observations. 
As for the Skyrme functionals, one suspects that they could not give the 
correct size of staggering as well. Such a conclusion may be supported by 
our DF3 calculations with parameters /| x = — f .79, = 1.66, = 0, i.e. 
with the pairing force of (4.5) which contains only a x 2 ^ 3 dependence, without 
gradient of density. The external strength constant /J x = —1.79 is taken from 
our preferable set (d), and the parameter = 1 .66 is found to get practically 
the same S n values as with the variant (d) itself. The results for the radii 
are shown in the lower panel of Fig. 6 by the dotted line. It is seen that the 
staggering appears but it is too small indeed. With a x 1//6 dependence, which is 
used with the SkP and SLy4 parametrization, the effect would be even smaller. 
Thus, to reproduce the scale of the odd-even effects in radii (and the scale of 
the kinks) within the local energy-density functional approach, the blocking 
mechanism should be enhanced by some peculiar density dependence of the 
effective pairing force, the linear dependence on p a with a=l , 1/3, 1/6... could 
not give the desirable size. 
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This fact clearly demonstrates the importance of including odd-mass nuclei in 
the fitting procedure of force constants. Now, in general, a complete description 
of the ground state of an odd nucleus is quite difficult which reflects itself e.g. 
in the difficulty to reproduce experimental magnetic moments. These depend 
sensitively on the time-reversal-odd components of the polarization induced by 
the odd particle. However, the operator r 2 is time-reversal-even, and the small, 
not-well-known, odd components of the polarization enter the expectation 
value only quadratically. If, with given interaction, (r 2 ) is not well described 
by a HFB-type ground state, there is no hope to reproduce it with a more 
sophisticated state vector. As can be seen from Figs. 5 and 6, it is possible 
to reproduce the order of magnitude of the staggering of (r 2 ) with density- 
dependent pairing force (4.3) containing a gradient term. 

Similar arguments apply to the nuclei off magic numbers. Looking only at even 
neutron isotopes, the slope of (r 2 ) as a function of iV changes at the magic 
numbers, producing a "kink". With simple interactions in the pp-channel, this 
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Fig. 7. Pairing gap A for 190 Pb (top) and 214 Pb (bottom) as a function of the radial 
coordinate r. The curves are marked in the same way as in Fig. 5. 



35 




Fig. 8. Pairing gap A for 190 Pb (top) and Pb (bottom) as a function of the 
dimensionless isoscalar density x. The curves are marked in the same way as in 
Fig. 5. 

is not reproduced in HFB-type calculations^. This relative increase of (r 2 ) off 
magic numbers is thought to be connected with "dynamical deformation" due 
to ground state fluctuations of the surface. These ground state fluctuations 
of the surface, being time- independent, actually can not be separated from a 
static diffuseness of the surface, and are present in an independent particle 
model too. In Section 6 below we will give some estimates of the contribution 
of RPA-type ground state correlations to this effect. It emerges that the main 
part of the effect should — and can — be obtained on the HFB level. 

In Fig. 7 the r-dependence of A for two isotopes 190 Pb and 214 Pb (the ends 
of the measured chain) is shown. One can see that with increasing external 

14 A kink has, however, been obtained in relativistic mean-field calculations, with 
simple pairing in the Pb chain [79] , at the expense of significantly too weak binding of 
the neutron single particle states above the Fermi level of 208 Pb which increases the 
polarization effects through neutron-proton interaction (see Ref. [80] for a detailed 
discussion of this point, and also our commentary to Figs. 9 and 10 of the present 
paper). 
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Fig. 9. Root-mean-squared neutron (n, open symbols) and proton (p, solid symbols) 
radii for lead isotopes calculated with functional DF3 and parameter set (d) of the 
pairing force (4.5) (circles), in comparison with the Skyrme-HFB calculations for 
the forces SLy4 [75] (squares) and SkP [74] (triangles), and with the RMF-HBCS 
calculations [81] (diamonds). The big open circle with error bars shows the "exper- 
imental" rms neutron radius in 208 Pb (see text). 

attraction oc /| x , A becomes smaller in the volume and gradually concentrates 
in the outer part of the nuclear surface. This is demonstrated more clearly 
in Fig. 8 were A is plotted as a function of dimensionless isoscalar density 
x. The curves (c)-(f) are crossing very nearly at the same point which has 
approximately the same coordinates in both cases: x piv ~ 0.7 and A piv pa 0.7- 
0.8 MeV. This "pivoting" point corresponds to the one seen in Fig. 1 for A in 
nuclear matter but in finite nuclei the value of A P i v turns out to be lower by 
a factor of 4 due to the repulsive gradient term oc /y. 

It is of interest to compare the predictions obtained with different mean-field 
approaches not only for the proton radii but also for the neutron ones. An 
example is given in Fig. 9 where the results for the lead isotopes obtained 
with the functional DF3(d) are shown in comparison with the HFB+SkP and 
HFB+SLy4 calculations, and with recently published relativistic mean-field 
Hartree calculations [81] based on the effective force NL3 and the constant 
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pairing gap prescription of Ref. [82] within the Bardeen-Cooper-Schriffer for- 
malism (RMF-HBCS). It is seen that the EDF model with gradient pairing 
yields staggering in the evolution of the proton and neutron radii, and pro- 
duces a kink in both of them at A=208. The Skyrme-type functional produce 
very small kinks in proton radii. The calculations with the RMF-HBCS model 
are given in Ref. [81] for even-even nuclei only, but one would guess that this 
model is not able to produce noticeable odd-even staggering too. At the same 
time, as seen in Fig. 9, both Skyrme-HFB and the RMF-HBCS functionals 
predict bigger neutron radii than our EDF calculations, and all of them yield 
a distinctive kink in the evolution of (r 2 ) n at A = 208. One observes rather 
close agreement between the proton radii obtained with all considered mod- 
els but a wide spread in the neutron radii. The RMF-HBCS neutron radii 
are particularly large. We remark at this point that the "experimental" rms 
neutron radius in 208 Pb of 5.593 fm deduced from the analysis of the high- 
energy polarized proton scattering in Ref. [83] is significantly lower than the 
prediction of the RMF(NL3) model; this fact has been already mentioned in 
Refs. [37,84]. But one should bear in mind that the result of this analysis 
is model-dependent, and the errors in the extracted neutron radii could be 
quite large. The difference between neutron and proton rms radii for 208 Pb, 
+0.14±0.04 fm, has been also deduced in Ref. [83], and we just put the cor- 
responding error bars for the (r 2 )}/ 2 value in this nucleus as shown in Fig. 9. 

Considering general trends in the behavior of a given sequence of nuclei, we 
notice that, because of the strong pn attraction, the evolution of the proton 
radius should be driven by that of the neutron one through self-consistency, 
and vice versa. Moreover, in neutron-rich nuclei, an anticorrelation should 
exists between neutron radii and neutron separation energies: the larger (r 2 ) n , 
the smaller S^n- One may suspect that too big neutron radii with a kink at 
N = 126, which is correlated with a kink in proton radii, would mean too 
small S 2n values beyond the shell closure. This is indeed the case as clearly 
demonstrated in Fig. 10. This figure displays two-neutron separation energies 
5*211 (upper panel) and differences of mean square charge radii S(r 2 h ) with 
respect to 208 Pb (lower panel) for the even lead isotopes. The results obtained 
with the RMF(NL3), Skyrme-SLy4 and EDF-DF3(d) models are compared 
there with experiment. It is seen that all these models reproduce values 
for lighter Pb isotopes below A = 208 fairly well, with more or less the same 
quality, the deviation from experiment being within 1 MeV or so (one should 
bear in mind that S^n for the A = 190-198 lead isotopes are derived in Ref. [70] 
from systematic trends, and the models discussed here describe these trends 
reasonable well). But for the heavier lead isotopes the predictions are different. 
Our EDF calculations with gradient pairing agree with experimental S 2n values 
within 0.5 MeV while Skyrme-SLy4 and RMF(NL3) models give deviations up 
to 1.5 and 2.0 MeV, respectively. It is also seen in the lower panel of Fig. 10 that 
the isotope shifts in charge radii for lighter Pb isotopes are nicely reproduced 
with our EDF approach, the 5{r 2 h ) values for the A > 208 nuclei and the kink 
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being slightly underestimated. The Skyrme-SLy4 and RMF(NL3) models work 
worse and yield results of equal quality for the lighter isotopes, but produce 
much different kinks and S{r^ h ) values beyond A = 208. The kink obtained 
with the RMF model incorporating simple BCS pairing, with constant pairing 
gap, is impressive indeed, and good description of this anomaly in the isotope 
shifts of Pb nuclei has been reputed as a considerable success of the relativistic 
mean-field approach. The above consideration, which has much in common 
with that of Ref. [80], points out, however, the importance of simultaneous 
description of both physical quantities, i.e. the energetic (SW) and geometrical 
( ^( r ch)) differential observables, in order to get a deeper insight into the nature 
of this anomaly in the isotope shifts. 




Fig. 10. Two-neutron separation energies (upper panel) and differences of mean 
squared charge radii with respect to 208 Pb (lower panel) for even lead isotopes. 
The open circles connected by the solid lines are obtained with functional DF3 and 
parameter set (d) of the pairing force (4.5). The open squares (diamonds) connected 
by the dotted (dashed) lines are from the Skyrme-SLy4 (RMF-HBCS) calculations. 
Solid circles: experimental data for S2n, including those derived from systematic 
trends, from [70], and data for S(r^ h ) from [71-73]. 
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Fig. 11. Neutron separation energies (top) and differences of mean square charge 
radii with respect to 116 Sn (bottom) for tin isotopes. Open circles: EDF calcu- 
lation with the set (d) of the pairing force, eq. (4.5). Open triangles (squares): 
Skyrme-HFB calculation with the force SkP (SLy4). In the lower panel, the open 
diamonds correspond to the RMF-HBCS S(r^ h }) values [81]. Solid circles: experi- 
mental data from [70] (for S n ) and [85] (for S(r^ h )). 

A=99 to A=136 are shown in Fig. 11 in comparison with experimental data 
and with the Skyrme-HFB calculations based on the SkP and SLy4 forces. In 
the lower panel of this figure, we also show the RMF-HBCS S(r^ h ) values [81]. 
The results for radii from these three models are available for even isotopes 
only. Of the six parameter sets (4.5) of the pairing force, only the results 
for our preferable set (d), with a scaling factor of 1.05, are shown. Similar 
to the lead chain, the S Q values for tin isotopes are reproduced reasonably 
well. However, on both ends of the considered chain, near the magic nuclei 
100 Sn and 132 Sn, the common trend for all mean-field calculations presented 
in the upper panel of Fig. 11 is that the size of the odd-even effect in neutron 
separation energies becomes noticeably smaller than the experimental one. As 
seen in the lower panel in Fig. 11, the evolution of the charge radii in the 
Skyrme-type and RMF-HBCS calculations is rather smooth, the SkP, SLy4 
and RMF functionals slightly overestimate <5(r 2 h ) in the heavier tin isotopes 
above A = 120. The desirable size of staggering in the behavior of <5(r 2 h ) for 
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Fig. 12. Two-neutron separation energies for even tin isotopes. Open circles: 
EDF-DF3 calculation with the set (d) of the pairing force, eq. (4.5), with a scaling 
factor of 1.05. Open triangles (squares): Skyrme-HFB calculation with the force SkP 
(SLy4). Open diamonds: RMF-HBCS calculation [81]. Solid circles: experimental 
data from [70]. 

tin isotopes could be obtained with gradient pairing as shown by open circles 
in this figure. At the same time, the EDF calculations with functional DF3 
significantly underestimate S(r^ h ) in lighter isotopes below A=114. This fact 
seems to be correlated with the weakening of the pairing when approaching 
A=100 and points out the need of both more careful adjustment of the normal 
isovector part of the energy-density functional and, probably, invoking the 
dependence of the pairing force on the isovector density p_ = p n — p p . 

The calculated two-neutron separation energies for even tin isotopes up to 
the neutron drip line are shown in Fig. 12. Again, the set (d) of the pair- 
ing force with a scaling factor of 1.05 was used in the EDF calculations with 
the parametrization DF3. One can observe noticeable differences between the 
Skyrme-HFB calculations with the SkP and SLy4 forces, the RMF-HBCS pre- 
dictions, and the EDF results. The RMF-HBCS calculations [81] are performed 
with deformed code, and available for the tin isotopes only up to A = 156; 
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Fig. 13. Isoscalar (top) and isovector (bottom) densities for some tin isotopes. 

these calculations produce some irregularities in the S 2jl values near A = 100 
and A = 132 (the nuclei beyond 140 Sn are predicted to be well deformed). At 
the same time, the three other spherical microscopic calculations give surpris- 
ingly similar predictions when approaching the neutron drip line, the heaviest 
bound tin nucleus being 172 Sn, 174 Sn and 176 Sn for the EDF, SkP and SLy4, 
respectively. One should notice that the Skyrme-HFB calculations were per- 
formed with a discretized continuum in a spherical box of finite size [74] which 
allows to artificially obtain a stable solution for the nuclear ground state even 
for positive value of the chemical potential. Such a solution is not possible 
with the coordinate-space technique used in the present paper because of the 
physical boundary conditions imposed for the scattering states. 

The evolution of the isoscalar and isovector densities along the tin isotope 
chain is illustrated in Fig. 13. It is seen that isovector density p_ = p n — p p 
increases both in the volume and at the surface as the neutron number N 
becomes larger. The isoscalar (matter) density tends to be slightly decreased 
with N in the volume and became more diffused at the surface. The matter 
distributions in the two magic nuclei, 100 Sn and 132 Sn, have a steeper slope 
at the surface compared to the non-magic isotopes. The radial dependence 
of the neutron pairing field A(r), of the anomalous density v{r) and of the 
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Fig. 14. Neutron normal density (top), anomalous density (middle) and pairing gap 
(bottom) in three selected tin isotopes. 

normal neutron density p n (r) for the stable nucleus 116 Sn, for the unstable 
isotope 150 Sn and for the drip-line nuclide 172 Sn are shown in Fig. 14. One 
can see that in the heavier isotopes, the anomalous density becomes more 
diffuse at the surface with a longer tail compared to the normal densities, 
in agreement with the discussion given in Appendix C. It is also seen that 
the pairing field A in the outer part of the nuclear surface has a prominent 
maximum which moves outwards and gets a larger tail near the drip line. The 
concentration of A outside of the half density point is related to the specific 
cancellation between the external attraction and the repulsive gradient term 
in the effective pairing force [24]. 
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Fig. 15. Neutron, proton and matter rms radii (top), mean energy of isovector 
dipole (middle) and Fermi charge-exchange (bottom) transitions for tin isotopes. 
The solid dots for (u;)di p correspond to the experimental data for the maximum of 
the dipole photoabsorption Lorentz curves, Refs. [86,87], while those for (cj)Fermi to 
the experimental positions of the isobaric analog states with respect to the daughter 
nuclei, Ref. [88]. 

The staggering phenomenon observed in the behavior of charge radii plotted 
as a function of neutron number is the one of the prominent odd-even effects 
which has been systematically measured and widely discussed. Similar effects 
are expected to exist for other quantities such as neutron and matter radii, 
centroid energies of multipole excitations (position of the giant resonances), 
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etc. 



Shown in the upper panel of Fig. 15 are the root-mean-square proton, neutron 
and matter radii for the tin isotope chain calculated with the functional DF3 
and with the set (d) of the pairing force. In the middle panel of this figure 
we show the mean energies (c^)dip of dipole isovector excitations. These are 
calculated as the square root of the ratio m^/mi with m\ and m 3 the linear 
and cubic energy-weighted sum rule, i.e. the first and the third moment of the 
corresponding RPA strength distribution, respectively [89]: 



(u) 



dip 



~3m NZ 



I 



dfdf ^JF"P(f,f') 9p P 



df 



dr' 



1/2 



(5.1) 



Here jF np is the effective neutron-proton interaction obtained from the energy- 
density functional as the second variational derivative S 2 Ei nt /Sp n Sp p . In the 
lower panel of Fig. 15 we show the mean energies (u) Fermi of the charge- 
exchange + excitations, i.e. the Fermi transitions, in tin nuclei with N > 
Z. These energies are calculated within the sum rule approach as the ratio 
(mf + mi )/(mQ — itlq ) with mX and m\ the first moment of the strength 
distribution of the Fermi transitions in the f3 + and /3~ channel, respectively 
(energy-weighted sum rules) and with itlq and itlq the corresponding non- 
energy weighted sum rules. The resulting expression for (u;) Fermi is given by 
(see, for example, [90]): 

(^)Fcrmi = N \ z J dr C/ C oui(r) {pjf) - p p (r)) , (5.2) 



where <7coui is the Coulomb mean field potential. 

It is seen in Fig. 15 that the rms neutron, proton and matter radii as functions 
of the mass number A reveal a kink at magic 132 Sn, and at larger A the dif- 
ference between rms neutron and proton radii starts to increase more rapidly. 
The staggering in radii is present mostly in the region between the two magic 
nuclei, from 100 Sn to 132 Sn, and this effect is practically washed out beyond 
A = 140. The mean energy of dipole transitions occurs to be in anticorrelation 
with such a behavior in radii, and this seem to be in qualitative agreement 
with experimental data (one should mention that eq. (5.1) overestimates the 
position of the giant dipole resonance by m 0.5 MeV because of the cubic 
energy-weighted sum rule m 3 used in the derivation of (cj) dip ). One also ob- 
serves a distinct kink in the behavior of (cj) dip at A = 132. Beyond this magic 
number the mean dipole energy decreases rather fast and then nearly saturates 
when approaching A = 172. This might be connected with an enhancement 
of the low-energy dipole transitions [91] and also with possible appearance of 
the so-called soft dipole mode [92] in nuclei near the neutron drip line. The 
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anticorrelations in the staggering behavior of (ou) dip and the rms radii {r 2 )l/ 2 , 
{r 2 )y 2 can easily be understood by considering the influence of pairing on 
the gradients of neutron and proton densities entering eq. (5.1). The odd-even 
effect in (cj)di P is constructive and more pronounced in the A < 132 region 
where both gradients strongly overlap. Beyond 132 Sn, the larger differences 
between neutron and proton rms radii imply lower overlap between neutron 
and proton density gradients at the nuclear surface, hence a smaller mean 
dipole energy. The situation with mean energy of the Fermi charge-exchange 
transitions is different. From eq. (5.2) one expects that the correlation between 
(cj) Fcr mi and the neutron and proton rms radii should be destructive. As seen 
in the lower panel in Fig. 15, the staggering in the evolution of (u;)Fermi with A 
is very weak indeed and almost invisible. The kink at the magic mass number 
A = 132 is also much less pronounced compared to the dipole case. Remark- 
able enough, the theoretical self-consistent sum-rule predictions for (a;) Fermi in 
tin isotopes are in excellent agreement with the available experimental data 
on the position of the isobaric analog states. 

Fig. 16 displays the neutron separation energies and isotope shifts in mean 
square charge radii for the calcium isotope chain. In this case, the EDF cal- 
culations with the preferable set (d), eq. (4.5), deduced from the Pb isotopes 
yield noticeable but rather small staggering both in S n and S(r 2 h ). A possi- 
ble way to get the desirable size of odd-even effects is to enhance the pairing 
correlations by introducing a scaling factor 1.35 for the pairing force (all pa- 
rameters being scaled in the same way). It is seen then in Fig. 16 that the 
neutron separation energies are reproduced fairly well, with more or less the 
same quality as with the Skyrme SkP and SLy4 functionals, though the effect 
of the scaling is rather mild. The effect on the behavior of charge radii is more 
drastic: the S{r 2 h ) values are increased by a factor of 3 to 4 yielding a very 
good agreement with experiment. This example demonstrates that the anoma- 
lous A-dependence of charge radii in Ca isotopes could be correctly described 
within the EDF approach with renormalized pairing force. However, the pos- 
sibility to determine a single parametrization of the effective interaction in the 
pp channel, which would be valid both for heavy and light nuclei, remains an 
open problem. In any case, we got an indication that the staggering in 5(r^ h ) 
of Ca isotopes, which is a challenge for nuclear models, could be explained 
with a peculiar density dependence of the pairing force. One can see in Fig. 16 
that the Skyrme SkP and SLy4 functionals produce a nearly smooth increase 
of radii with A, with weak irregularities below and above A = 48, but still 
with rather sizeable kinks at double magic 48 Ca. The RMF-HBCS predictions 
for 5(r 2 h )i Ref. [81], are much different: the charge radius is nearly a constant 
up to 48 Ca, with a sudden increase for 50 Ca. Interestingly, the functional DF3 
yields, in agreement with experiment, a small negative isotope shift in charge 
radii between 40 Ca and 48 Ca, which does not depend on pairing. 

The EDF calculations, as demonstrated in Fig. 17, show that the trends and 
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Fig. 16. Neutron separation energies (top) and differences of mean squared charge 
radii with respect to 40 Ca (bottom) for calcium isotopes. Open circles: the EDF 
calculation with the set (d) of the pairing force, eq. (4.5). Solid lines connect the 
EDF results also obtained for the set (d) but with a scaling factor of 1.35. Open 
triangles (squares): Skyrme-HFB calculation with the force SkP (SLy4). In the lower 
panel, the open diamonds correspond to the RMF-HBCS 5(r^ h ) values [81] (for even 
isotopes only). Solid circles: experimental data from [70] (for S n ) and [93-95] (for 

the staggering effects in the mean energy of isovector dipole excitations and 
charge-exchange Fermi transitions for calcium isotopes are in strong anticor- 
relation with evolution of the ground state mean squared radii. Similar to 
the case of the tin isotopes discussed above, the anticorrelation in (c^Fermi 
is weaker than in (a;) dip- Unfortunately, there are no systematic experimental 
data on dipole excitations in Ca isotopes at our disposal, from which the (a>)di p 
values could be deduced. As for the mean energies of the Fermi transitions, 
which are calculated with the self-consistent EDF sum rule approach and plot- 
ted as a function of A in the lower panel in Fig. 17, they may be compared 
with the experimental positions of the isobaric analog states. It is seen that 
the sum rule expression, eq. (5.2), underestimates the energy of the analog 
states in calcium isotopes by some 200 keV but reproduces qualitatively the 
observed trend. 
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Fig. 17. Neutron, proton and matter radii (top), mean energy of isovector dipole 
(middle) and Fermi charge-exchange (bottom) transitions for calcium isotopes. The 
solid dots for (w)Fermi correspond to the experimental positions of the isobaric analog 
states with respect to the daughter nuclei, Ref. [88]. 



As, in the past, RPA ground state correlations have been invoked to explain 
details of the isotopic dependence of nuclear charge radii, we estimate the 
influence of RPA ground state correlations. 
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6 Ground state fluctuations 



An important point is the issue of RPA-type ground state correlations. In 
several papers [72,96,97] the influence of the average mean square deformations 
on the mean square charge radii has been discussed in a phenomenological way 
which is based on the following sum rule: Take the identity for the multipole 
operator = J dfr x Yx^(fl)ip^ (r)ip(r) 

E(0MQa,/|0) = EIWaJ0)| 2 . (6.1) 

The right hand side is equal to the summed transition strength X) s -B(EA;0— >s), 
the left hand side is the expectation value of the squared multipole operator 
Qx in the ground state. In semiclassical nuclear models, this expectation value 
is proportional to the square of the deformation of multipolarity A. Thus, even 
for spherical nuclei the rms deformation does not vanish. In the liquid drop 
model this is interpreted as due to zero point surface vibrations, but it is to be 
emphasized that there is no observable time dependence, and no breaking of 
the spherical symmetry due to this "dynamical" deformation. Experimentally, 
the effect can not be distinguished from diffuseness of the surface arising from 
any other mechanism. 

Our method (as well as any HF or HFB calculation) incorporates only the dy- 
namical deformation effects corresponding to a gas of noninteracting particles 
(magic nuclei) or of a gas with only pairing (nonmagic nuclei). To include the 
RPA ground state correlations, the calculated charge radii should be increased 
by 

Hrl) = 5(r c 2 h ) RPA - 5(r 2 ch ) UFB . (6.2) 

However, the parameters have been adjusted to reproduce the experimental 
density distribution and therefore the average RPA correlations are included 
due to the choice of parameters. In a high precision fit, this should be corrected 
for, by subtracting the value given by (6.2) for the reference nucleus. 

A microscopic method to calculate the corrections to (r^ h ) due to low-energy 
surface vibrations has been proposed in [98]. We use the liquid drop model 
and eq. (6.1) to estimate the influence of isoscalar collective excitations on the 
charge radii, as proposed by Esbensen and Bertsch [96]. 

As already noted in [96], the inclusion of noncollective states in eq. (6.1), which 
can not be considered to be surface oscillations will not do much harm, because 
the contribution of these states will be essentially the same to 5{r^ h ) RPA and 
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^( r ch)HFB and will drop out. According to [72,99,100], for multipolarity A > 

Hrl) ee (r 2 ) d - (r 2 ) s = ^<r 2 > s £<^> , (6-3) 

A 

where the suffixes d and s mean (dynamically) deformed and spherical, re- 
spectively. 

Within the same model 

S°(EA) ee ^2 B(EX; O^s) = £<0| \Q X ^\0) = (^^f <#> ■ (6.4) 

With the help of (6.4) and (6.1) the sum over (3\ in eq. (6.3) can be replaced by 
the sum over the total B(EX) -values, and we reach the following expression: 

4-ir 1 

Hrl) = E J^S [^rpa(EA) - 5° FB (EA)] . (6.5) 

A 

Of course, for A = 1 one has to exclude the spurious translation. For A = we 
have to use a compressible drop model, and obtain similarly 

^ = 108zSiZg [g ^ A(E0) ~ ^fb(EO)] . (6.6) 

where S'(EO) is calculated according to eq. (6.4) with Q 00 = r 2 F 00 . 

The excited states have been calculated in quasiparticle-RPA (QRPA) which 
has been used in the form (we adopt the same symbolic notations as in 
Ref. [10]): 

V = e q V + FAV , S° PA = --lm(e q V Av) , (6.7) 

7T 

where e q is the quasiparticle local charge with respect to the external field 
Vo- For a long-range electric field, Vq oc Q\ tlJl , from gauge invariance one gets 
e q — 1 [10]. The A matrix is written in detail in [101]. 

As in our formalism the effective interaction amplitudes F are obtained as 
the second functional derivatives of E- mt with respect to the corresponding 
densities, one can see here that due to p-dependence of £ a nomai, eq. (4.1), there 
are mixed terms for the effective interaction, F^ and F&, which do not vanish, 
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Table 1. Excitation energies (MeV) and transition probabilities (e 2 •fm 2i ) for low- 
lying collective states. Experimental data are taken from Refs. [102,103] (for Ca), 
[104] (for Sn) and [105] (for Pb). 
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and which couple the ph with the pp or hh channel (r-indices are neglected 
for simplicity): 

= 62E ;f U] e—i([ P , u]; f) df , (6.8) 

opov opov J 

A comparison of QRPA results for the first collective 2 + and 3 _ states with 
experiment is presented in Table 1 for selected nuclei from the Ca, Sn and Pb 
isotope chains. The influence of the spin-orbit part in the ph interaction has 
been investigated in [21]; the results given here have been obtained without 
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Table 2. S(r^ h ) (in fm 2 ) for different multipolarities for 40 Ca, 44 Ca, and 48 Ca iso- 
topes calculated within the framework of self consistent QRPA. 
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this spin-orbit contribution. Comparison with [21] shows that here the changes 
in the functional did not improve the agreement with experimental data. 

We restricted our consideration to 2 + and 3 _ states for all isotope chains 
because the other multipolarities contribute much less to the mean square 
charge radii [98]. This has been checked by computing also S(r^ h ) from + and 
4 + states for three calcium isotopes. In all these calculations the functional 
parameter set of Ref. [22] have been used. These results are presented in 
Table 2 in comparison with the calculations for 2 + and 3~ multipolarities. One 
can see a here significant difference between 2 + and 3 _ states on one hand and 
+ and 4 + states on the other hand. Therefore it is sufficient to include only 
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Fig. 18. Contribution of RPA ground state correlations to the squared charge radii in 
calcium isotopes. Shown are the 2 + (diamonds) and the 3~ (triangles) contributions 
and the sum of both (circles). It is seen that 40 Ca is quite soft for 3~ deformation 
and gets a big contribution from the corresponding ground state fluctuations. 
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Fig. 20. The same as Fig. 19 for lead. 

strong phonons which contribute significantly to the mean square radii. This 
result is very similar to that of Barranco and Broglia [97] which they obtained 
with the approach of Esbensen and Bertsch [96] for the calculation of S{r^ h ) 
values. However, it should be mentioned that they used a non self consistent 
model with more freedom in the choice of parameters. 

Let us now look at the results of our calculation for the calcium, tin and lead 
isotope chains which are presented in Figs. 18, 19, and 20, respectively. One 
can see here that for calcium isotopes 3~ states give the main contribution to 
^( r ch)tot, especially for 40 Ca (Fig. 18). The 2+ contribution is not very big, in 
contrast to the result of [97]. This is connected with the fact that we could 
not obtain the experimental values for excitation energies and especially for 
B(E2) values in calcium isotopes in the middle of the chain. Presumably for 
good reproduction of the experimental values one has to go beyond the QRPA 
in a model which allows to take into account more complex configurations 
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than lplh. Nevertheless, it is shown that there is considerable influence of 
ground state correlations on the calcium mean square charge radii, at least in 
second order in the phonon amplitudes. Combining the EDF results shown in 
Fig. 18 with those in Fig. 16 one sees that the staggering in charge radii of 
calcium isotopes could be enhanced by the calculated ground state correlation 
corrections; however, this would produce too big negative isotope shift between 
40 Ca and 48 Ca thus deteriorating the good agreement with experiment seen in 
Fig. 16 alone. 

In Figs. 19 and 20 the ground state correlation contribution to S(r^ h ) is shown 
for a few isotopes of the tin and lead chains. These nuclei are more stiff, 
and QRPA works much better (see Table 1). The results show that here the 
ground state correlation contribution is only a small fraction of the value 
obtained without correlations. Moreover, the general trends look roughly like 
the trends of the curves for 5(r^ h ) from Figs. 5, 6 and 11, which therefore are 
enhanced but not qualitatively changed by the correlations. 



7 Summary and conclusions 

In this work the energy density functional method for superfluid nuclei has 
been described in detail. The generalized variational principle has been formu- 
lated for the local functional with a contact effective force T^- in the particle- 
particle channel and cutoff energy e c > e-p. For spherical finite systems, the 
coordinate-space technique, involving the integration in the complex energy 
plane of the Green's functions obtained by solving the Gor'kov equations ex- 
actly, with physical boundary conditions both for bound and scattering quasi- 
particle states has been used. The technique has been extended to odd nuclei 
by using a uniform filling approximation. 

Within the framework of this EDF approach, the combined analysis of dif- 
ferential observables including odd-even mass differences and odd-even effects 
in charge radii was performed for a few isotopic chains. The staggering and 
kinks in 5(r^ h ) were shown to be very sensitive to the density dependence of 
the effective pairing force. The observed isotopic shifts were successfully de- 
scribed with the force containing a density-gradient pairing term. With such 
a parametrization, as the external attraction oc /| x increases, the pairing field 
A(r) becomes smaller in the volume and concentrates in the outer part of the 
nuclear surface. 

A few possible parameter sets extracted from the lead chain were used to study 
the ground state properties of uniform nuclear matter with s-wave pairing, in 
particular the behavior of the energy gap at the Fermi level as a function of 
density. The deduced sets, which give a satisfactory description of both the 
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neutron separation energies S Q and isotope shifts S(r^ h ) in charge radii, yielded 
about the same value of A pa 3.3 MeV at k-p = 1.15 fm _1 (at pa 0.65 of the 
equilibrium density). It was found that some sets in the dilute limit correspond 
to the strong coupling regime. The situation when the regime changes from 
weak to strong pairing was considered in detail, and, in strong coupling, the 
properties of dilute matter were shown to be completely determined in leading 
order by the singlet scattering length a nn . In weak coupling, near gap closure, 
the pairing gap Ap is shown to be directly expressed through the Fermi level 
phase shifts. The density-dependent cutoff contact pairing interaction corre- 
sponding to the realistic value of the singlet NN scattering length was found 
to be a preferable choice between the deduced parameter sets to be used in 
the EDF calculations. 

The isotopic chain of tin isotopes from the proton to the neutron drip line was 
calculated, and the evolution of the normal and anomalous (pairing) density 
distributions with A was analyzed. The parametrization DF3 of the normal 
part of the functional and the density- dependent pairing interaction provided 
fairly good description of odd-even effects in S n and S{r^ h ); however, in the 
region of lighter tin isotopes the slope in 5{r^ h ) plotted as a function of the 
mass number A come out steeper than in experiment. This might indicate 
that dependence on isovector density p_ = p n — p p should be incorporated 
in the effective pairing interaction. Anomalous behavior of charge radii in Ca 
isotopes was reproduced in excellent agreement with experiment within the 
EDF approach by scaling the pairing effective interaction which was extracted 
from the lead chain by a factor of 1.35. 

It was shown with the self-consistent sum rule approach that the density- 
dependent pairing induces sizeable staggering and kinks in the evolution of 
the mean energies of multipole excitations along isotopic chains. These effects 
were found to be in anticorrelation with the behavior of the ground state radii. 

Using quasiparticle-RPA multipole strength distributions obtained with the 
self-consistent interaction, the contribution of phonons (ground state correla- 
tions) to the differences of charge radii was calculated by the method proposed 
by Esbensen and Bertsch [96]. It was shown that these phonon corrections, at 
least for tin and lead isotopes, are quite small in comparison with the values 
provided by the main HFB contribution, and do not lead to qualitatively new 
effects. 

The results obtained in the present paper were compared with the conventional 
Skyrme-HFB, Gogny-HFB and relativistic mean field Hartree-BCS calcula- 
tions. On the whole, a better description of the differential observables was 
achieved with the EDF method, the most important effect which was suc- 
cessfully reproduced is the staggering in charge radii. Up to now, this effect 
was practically missed in all other mean-field models. The famous kink in the 
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isotope shifts of Pb isotopes has been also well described by our EDF calcu- 
lations, with density-dependent gradient pairing, and the crucial point is that 
the neutron separation energies are reproduced as well. These results support 
the conclusion drawn by Reinhard and Flocard [80] that this is apparently 
not the case in the RMF model. These authors showed that the description 
of the kink in Pb chain is not proprietary to the relativistic model, it can 
be also produced with the Skyrme functionals by generalizing the spin-orbit 
term. However, even with their "best" set SkI4, Reinhard and Flocard have 
found a large gap in the single-particle neutron spectrum similar to that of 
the RMF model. Our conclusion is that the physics behind the anomaly in 
the Pb radii is not solely due to the spin-orbit (which, moreover, has nothing 
to do with the staggering!), and that the simultaneous reproduction of both 
the energetic and geometrical observables is important to settle the problem. 

Unfortunately, a universal parametrization of the pairing force is still lacking, 
and for various isotope chains, especially for lighter nuclei, the parameters of 
had to be chosen somewhat different. It seems to be necessary to use a more 
refined parametrization, e.g. with additional dependence of on p_. Also, 
more attention should be payed to the normal part of the density functional 
in order to choose more carefully its parameters and improve the description 
of the bulk properties of nuclei such as their masses and absolute values of 
radii [106,107]. These questions will be investigated in future work. 
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Appendix A. General variational principle and ground state energy 
in local approximation 



In this appendix we show that, in the case of weak pairing |A| <C ep, the 
general variational principle for a local density functional with an energy cutoff 
e c > eF can be used to describe the ground state properties of superfluid 
systems. This allows to perform a full HFB-like calculation with a renormalized 
gap equation and an effective pairing 5-force. 

We start with the usual definitions and equations which follow from the general 
variational principle. The total energy E of a superfluid system, as given by 
eqs. (2.10)-(2.12), is a functional of the generalized one-body density matrix 
R. This matrix can be expressed in terms of the ground state expectation 
values of the pair products of the Landau-Migdal quasiparticle creation and 
annihilation operators ip: 



£(1,2) 



/ <^(2Mi)> (mm) 

\ <^W(l)> (V(2)V f (l)> 
p(l,2) v{\,2) 

2) <5(l,2)-p*(l,2) 



(A.l) 



Here the coordinate space representation is used. The numbers in brackets 
stand for the set of all relevant arguments and indices, e.g. (1) = ({r, s z , r 3 }i). 
In the following, / dl means integration over the continuous and summation 
over the discrete variables of this set. 

The Bogolyubov quasiparticle creation and annihilation operators fy, f3 a are 
defined by the transformation 




u* a (i) vm \ i m 

V a (l) U a (l) Lt(i) 



(A.2) 



This can be inverted to express the Landau-Migdal quasiparticle operators in 
terms of those of Bogolyubov 

« U„(i) u-(i) I, at 



The generalized Bogolyubov quasiparticle density matrix Q is given by 

Qaa ' ~ I ) {pM ) ■ {AA) 
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By demanding now the ground sate |$o) to be the Bogolyubov quasiparticle 
vacuum, (3 a \®o) — 0, the matrix Q becomes diagonal in the a-representation: 

Qaa'=8aa>Q, with Q = I I , (A. 5) 

1 



and for (A.l) one obtains a "supermatrix" formula: 

3(1,2) = EW a (l)Q m) }Vi,(2) =£W a (l)QMi(2) , (A.6) 



with the matrix W of the Bogolyubov transformation (A. 3): 

W Q (1) 



VK(1) C/*(l) 



For the generalized density matrix we thus find 



(A.7) 



6/ 19) _r/W) J£(l)0a(2)\ , , 

l ' J ^U*(1)K(2) U* a (l)U a (2) I ■ 



The matrix W obeys the generalized closure and orthogonality relations: 

E W«(l)Vl£(2) = 1 ° W 2) ee 75(1, 2) , (A.9) 
^0 1 J 

J dlWt(l)>V Q ,(l)=/(J aa ,. (A.10) 

The generalized density matrix R is Hermitian and from eqs. (A.5-A.10) fol- 
lows its important property (2.15) which, in the coordinate space representa- 
tion, reads 

f d3#(l, 3)3(3, 2) = R(l, 2) . (A.ll) 



This property reflects the fact that the ground state is a HFB quasiparticle 
vacuum. 

Minimizing the energy and imposing the average particle number conservation 
(2.14) and the relation (A.ll) as constraints, we arrive at the variational prin- 
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ciple 5I[R] 
here as 



= with the variational functional of eq. (2.16) which we rewrite 
I[R] = E[R] - /iTrp - Tr(A(_R - R 2 )) . (A.12) 



The trace of a supermatrix A, in the coordinate space representation, is defined 
by 



, A 11 A 12 1 r ^ .. 
Tr^Tr( i2ii22 ]^/dlE^(M). (A.13) 



The total energy of the system, eqs. (2.10)-(2.12), is the functional 

E[R] = Tr(tp) + E int [R] = Tr(tp) + S int(normal) [p] + E anomal [R] 

= E noimal [p] + E anoma i[R] , (A. 14) 

where t is the free kinetic energy operator and E int [R] is the interaction energy 
containing both normal and anomalous terms. Because of condition (A. 10), 
the matrix A act i, in the Bogolyubov quasiparticle space, may be taken to be 
diagonal: 



A. aa i 



1 ' E a " ' 
V E a j 



5 aa , = IE a 5 aa , , (A. 15) 



where E a is the set of Lagrange multipliers. In coordinate space representation 
one gets 

A(l,2)=5>aW a (l)Vl2(2). (A .16) 



The variational equation reads 

5 (E\p, v\ - /iTrp - Tr(A(_R - R 2 ))) = . (A. 17) 

For the variation of the first two terms in this equation we have 

5(E[p, v] - pVip) = Tr(H5R) = J dl d2W(l, 2)5^(2, 1) , (A. 18) 
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where the effective quasiparticle Hamiltonian Ti is given by 

^ / Ml, 2) -Ml, 2) A(l,2) \ 

n[i > Z) ~ { -A*(l, 2) 2) - h*(l, 2) J ' (A - i9j 



with 



Ml 21-^1 A(l 2 )- ^" ] - 5EM (A20) 



The variation of the generalized density matrix, 



572(1,2) = ^ ^ , (A.21) 

-<&/•(!, 2) 2) 



may be written in the form 

5/2(1,2) =E W «( 1 )^-' W 1'(2), (A.22) 



where SQ aa > is the variation of the Bogolyubov quasiparticle density ma- 
trix (A. 4). The choice of the matrix elements 5Q aa i as independent variables 
leads to the HFB equations in a simple way. The variation of (A.ll) gives 



Eq.(A.18) may thus be written as 



'i 



-1 



SQaa>Wi(2). (A.23) 



act' 

= £/ dld2^WtJ(l)^ c (l,2)W^(2)5Q^, (A.24) 

act' abed 

and the variation of the last term in (A. 17) may be expressed as 

E \ 

5TrA(R-R 2 )) = £Tr I a \SQaa. (A.25) 

-E a 
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Thus, the variation equation (A. 17) becomes 



ETr 



yvlnw a - 5 a , a 



E a o 
-E a 



8Qaa' — 



(A.26) 



Since the elements 5Q aa r are independent, we obtain the matrix equation 



E a 
-E n 



(A.27) 



which is equivalent to the HFB equations 
d2W(l,2) 



J d2W(l,2) 



[ U a (2) 


)-«.( 


Ua(l)\ 






V a (l) ) 


V a *(2)\ 


= — E a 


V Q *(i) N 









(A.28) 



(A.29) 



The matrix of second variational derivatives of the energy functional is the 
supermatrix of the effective quasiparticle interaction amplitudes T . In the 
pp-channel one finds 



JF pp (l,2;3,4) 



*A(1,2) 
5i/(3,4) 



(A.30) 



We assume that JF PP is a functional of the normal quasiparticle density ma- 
trix p and that the term -E an omai m (A. 14), which depends explicitly on the 
anomalous density matrix, may be written as 

E„4R] = \ J dl-d4z/t( 2 ,l)^PP(l,2;3,4;[p])z/(3,4), (A.31) 



with the antisymmetrized pp-interaction 

^(1, 2; 3, 4) = 2; 3, 4) - 2; 4, 3) 

The gap equation reads 

A(l,2) = || d3d4^(l,2;3,4;[p]M3,4). 



(A.32) 



(A.33) 
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The expression for the anomalous energy (A. 31) may thus be written in the 
form 

^anomal = \ J dl d2z/ t (2, 1)A(1, 2) . (A.34) 

To proceed further and perform the renormalization procedure it is convenient 
to use the Green's function formalism (see Appendix B). From (A. 8) and 
eqs. (B.9)-(B.ll) it follows that the generalized density matrix R may be 
represented by a contour integral of the generalized Green's function G in the 
complex energy plane: 

= / c ^G(l,2;e). (A.35) 



The integral is performed along a contour C which goes from — oo to below 
the real e-axis (Ime = —0), crosses this axis at e = and goes back to — oo 
above it (Ime = +0), with the energy variable e measured from the chemical 
potential //. The contour is shown schematically in Fig. 21, see Appendix C. 

The normal and anomalous density matrices are then given by 

P(l,2)= / ^-G s (l,2;e) = £ ^(1)^(2) , (A.36) 

JC 2711 (E a >0) 

"(U) = f ^F(l,2;e)= £ l£(l)tf a (2) . (A.37) 



The summation (or integration) in these expressions extends, in principle, to 
infinity, the convergence being dependent on the properties of the effective 
interaction. It is well known that in the case of a local force the anomalous 
density matrix z/(l, 2) diverges at r\ — > r 2 . In the gap equation (A. 33) with v 
expressed by (A.37), the integration over the region far from the Fermi sur- 
face can be avoided by using the standard renormalization procedure with 
an arbitrary energy cutoff [10]. But to calculate the energy of the system, in 
particular its (negative) anomalous part (A.34), one needs to know the "un- 
truncated" density matrices. On the other hand, as it is also well known, the 
total pairing energy E pair contains a (positive) contribution from the kinetic 
energy term so that the sum of the two contributions converges more rapidly 
and -Ep a i r gets a finite value even in the case of a local pairing field (at least 
in infinite matter, see Appendix B). We shall show that the energy density of 
the uniform infinite system may be calculated exactly in the leading order in 
A 2 (p F )/ep by using a renormalized gap equation with an arbitrary cutoff e c 
but such that e c > ep, and by applying the variational principle SI[R] = to 
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the functional I[R] given by (A. 12) but considered now as a functional of the 
corresponding cutoff generalized density matrix R c . 

We introduce an energy cutoff e c > e-p. The contour C is split into two parts 
Ic = Ic c + Ic B i where C c lies below and above the real e-axis for — e c < Re e < 
and the contour Cb lies below and above it for Ree < — e c . Then we have 

z/(l,2) = z/ c (l,2) + 5 c z/(l,2), (A.38) 



where 



1,2)= J2 K*(l)f4(2)= / ^F(l,2;e), (A.39) 



0<E a <e 



and 



The renormalized gap equation may be written in the form 

A(l,2) = | J d3d4Ff(l,2;3,4K(3,4). (A.41) 

where the effective scattering amplitude is to be found from the equation 

^(1,2; 3, 4) = ^(1,2; 3, 4) 

+ J d5---d8J !rPP (l,2;5,6)fi(5,6;7,8)J !rf (7,8;3,4), (A.42) 

with 

B(l, 2; 3, 4) = / p:G (l, 3; e)G a (4, 2; e) . (A.43) 
JC B Alii 

Thus the anomalous energy (A. 34) may be represented by a sum of two terms: 

-E'anomal -^anomal ^c-Eanomal 

+ ll ili2 IcM F,{2 ^ )Ail - 2) - (A ' 44) 
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The normal density matrix is written in a similar way as a sum 

p(l,2)=p c (l,2) + Ml,2), (A.45) 



where 

i,; (1)1, ,(2) = / <l< 

0<E a <e c 



5 (1,2)= E KV)V a (2)=l ^G.(l,2;c), (A.46) 



and 



Ml, 2)= £ l£(l)V«(2)= / A Gs( i )2;e ). (A.47) 



E a >e, 



Using the definition of the single particle Hamiltonian h, eq. (A. 20), the 
amount of the total energy related to 5 c p in first order may be found by 
varying the normal part of the energy functional with respect to the normal 
density: 

(^normal = J dl d2/l(l, 2)S cP (2, 1) . (A.48) 

This includes also, if T vv depends on p, a contribution from the variation of 

-E'anomal • 

With 5 c -E a nomai from (A. 44) and S c E noima i defined by (A.48) we obtain the total 
change of the variational functional related to the cutoff: 

5 C (E - pN) = J ^-XhG s + Va - fj,G s ) . (A.49) 

This can be written in another form using the relation 

AF f = (e - h + n)G s - 1 (A.50) 



which follows from the Gor'kov equations (B.3): 



6. (E - iiNtji)) = ilY J |i( e + h - »)G. 

To get the second equality in this equation we have used the relation 

Tr J c de(e + h — p)G = which holds because the Green's function G 
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is diagonal in the basis of eigenfunctions of h (see eq. (B.8)) and has no 
singularities embraced by the contour Cb- 

In lowest order in A, from (B.3) one has = G A*G and we get 



5 C (E - /iiV(/i)) = ±Tr [ ^(e + h- fi)G AG A*G 

Z JCr iTll 



(cj - e fc )|A 



|2 



_ «5*U5*(«.'+*-i)V" (A,52) 

Here A^ are the matrix elements of A in the basis of the single-particle Hamil- 
tonian h with being its eigenvalues (see Appendix B, eq. (B.7)). 

In uniform infinite matter, the pairing field A(l, 2) is a function of |fi — r 2 |, all 
its nondiagonal matrix elements vanish and we see that 5 C (E — fj,N(/i)) = in 
the first order upon imposing the energy cutoff e c > cf- In finite systems, the 
change of S c (E — fiN(fi)) vanishes in the first order in diagonal approximation. 
A nonvanishing right-hand side of eq. (A.52) arises only due to nonuniformity 
of A, which means that in finite systems with more or less uniform density in 
the volume (like heavy atomic nuclei), this would be a surface effect. Reliable 
estimates of the sum in (A.52) involving only nondiagonal matrix elements 
Ajfc i i-^k are not easy to obtain. They depend strongly on the actual distribution 
of the pairing field in real nuclei, in particular in the surface region. In larger 
systems, with presence of volume pairing, the nondiagonal matrix elements of 
A are relatively small (~ A -1 / 3 , or less) so that the change of E — /iN due to 
the cutoff is expected to be at least by a factor of A^ 1 ^ 3 smaller compared to 
-Ep a i r , and therefore we can write 

\5 C E - fi5 c N(n)\ < |£ pair |, (A.53) 



where 

5 c N(fi) = J dlMl, 1) = N(ijl) - NM , (A.54) 

with N c the number of particles corresponding to the cutoff density p c of 
eq. (A.46): 

N c = J dlp c (l,l). (A.55) 

The pairing energy E pa i r is given in Appendix B by an exact expression (B.41) 
and by a BCS estimate (B.46). 
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The change in the system energy related with the cutoff may be written as 

5 C E — E[R] — E[R C ] = p5 c N , (A. 56) 



with S C N defined in (A. 54). Minimizing the energy functional with a fixed 
cutoff e c would give an additional change 

5E C = E[R C ] - E C [R C ] = p c 5N c , (A.57) 

where SN C is the change in the particle number within the cutoff energy space 
e < e c . In first order we have 

E[R\ = E C [R C ] + p5 c N + p c 5N c = E C [R C ] + p(5 c N + 8N C ) . (A.58) 

Of course, the constraint Trp c = has also to be used for the cutoff functional 
which gives 5 C N + 5N C = 0. 

We thus come to the important conclusion that passing to the cutoff functional 
with e c > ep leaves, to the first order in the weak pairing approximation, the 
energy, the variational functional E — fiN and the chemical potential /i of the 
system unchanged (to the extent that the estimate (A. 53) for a given nucleus 
is correct). 

According to the Hohenberg-Kohn theorem [27], the nuclear ground state 
properties can be described by a local density functional. One may there- 
fore assume that the anomalous energy is a functional of the local density p 
too. The above consideration implies that this should be also valid for the 
cutoff functional. Both the cutoff anomalous density and the effective pairing 
interaction JF^ may then be regarded as local functionals of p. We introduce 
a local approximation to JF^: 

^(l,2;3,4)=^(f 1 ,r 1 ) 

x S(n - f 2 )<J(fi - f 3 )5(f 2 - r A )8 SlSz 8 S2SA 8 TlTz 8 T2T4 . (A.59) 

where 

J*(?, T ) = ^([p c (f,n),p c (f,p)];r) (A.60) 

is a local functional of the cutoff quasiparticle density 

p c (f 1 ,T 1 )=Tr ai J d2p c (l, 2)5(2,1), n G {n,p} , (A.61) 

where 5(1,2) = 8(n - f 2 )8 SlS2 5 TlT2 . 
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In the case of a local functional it is convenient to work with the so-called 
anomalous density in the real space. For the time-reversed single-particle states 
\k) = T\k), in the coordinate space representation, we define 

&(1) = J d2T(l, 2)0*(2) , 0|(1) = J d20 fc (2)Tt(2, 1) , (A.62) 
with the operator 

T(l,2) = 5(f 1 -r 2 )(-l) 1 ^5- SlS2 5 TlT2 . (A.63) 

This operator is antisymmetric, T(2, 1) = — T(l,2), and has the properties 
/ d3T(l, 3)T(3, 2) = -5(1, 2) and J d3T(l, 3)Tt(3, 2) = 5(1, 2). The anoma- 
lous density is then defined by 

u c (f 1 ,T 1 ) = iTr ai J d2z/ c (l,2)Tt(2,l), (A.64) 

and its complex conjugate by 

v* c (n, n) = ±Tr sl | d2T(l, 2)^(2, 1) . (A.65) 

In the local approximation we get 

A(l,2) = T(l,2)A(f 2 ,r 2 ). (A.66) 

With these definitions, the anomalous part of the cutoff functional reduces to 
the expression (2.25), and for the local pairing field A(r, r) one obtains the 
gap equation in the simple multiplicative form (2.27). 

One may notice that our approach does not imply a cutoff of the basis since the 
general variational principle can be formulated with a "cutoff" local-density 
functional from which the ground state characteristics of a superfluid system 
may be calculated through the solutions of the Bogolyubov equations (A. 28) 
at the stationary point. To construct the normal and anomalous densities, en- 
tering this local functional, only those solutions from the whole set are needed 
which correspond to the eigenenergies E a of the HFB Hamiltonian up to the 
cutofip 5 ] e c > ep. This means that one can implement a local pairing effec- 
tive force in the HFB (or Gor'kov) equations and, most important, in the gap 

15 At first glance the situation looks similar to the HF case with a EDF without 
pairing where only the occupied orbitals up to ep in the self-consistent mean field are 
needed to calculate the density (a "natural" HF cutoff). But it should be emphasized 
that the Bogolyubov solutions cannot be obtained from the HF ones by perturbative 
methods, they spread to infinity in the energy space at any nonvanishing gap, hence 
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equation, renormalized with the same energy cutoff e c , in which the anomalous 
density by construction does not diverge. The whole set of resulting equations 
can be solved, for spherical nuclei, exactly with a coordinate-space technique 
(see Appendix C). It should nevertheless be stressed that the main difficulty 
now is to find such a local effective density-dependent pairing interaction 
that would be an universal one and would allow to reproduce the pairing 
ground-state nuclear properties (at least for a large domain of heavy nuclei). 



Appendix B. Gor'kov equations and pairing energy 



In this appendix we derive the general expression for the pairing energy using 
the Green's function formalism. We show that, in the case of a weak local 
pairing field, in infinite matter this expression reduces to the BCS formula 
used in eq. (4.30). 

We start with Gor'kov equations [44] for the generalized Green's function G(e) 
written in matrix form, 

(e-H)G(e) = I, (B.l) 



where H is the effective single-quasiparticle Hamiltonian (A. 19) and 

(G,(t) F(t) 



G < £ >=^(e) <B ' 2) 



with G s the normal and F the anomalous Green's function, respectively. For 
G s and F\ eq. (B.l) yields the set of two equations 

G s (e) = G (e) + G (e)AFt(e) , F\e) = -G (e)A*G s (e) , (B.3) 



where G and G are the single particle Green's functions connected with the 
normal parts h — /i and /i — h* of the effective Hamiltonian, respectively: 

Go(e) = i— , G (e) = —1 . (B.4) 

e — h + (/, e + h* — ji 



the convergence problem for the densities. This problem, as usually thought of, can 
be avoided only if one uses a finite-range effective force ensuring the convergence of 
the anomalous density at some energy which depends on the range of nonlocality of 
the force. Then the HFB equations become nonlocal leading to a technical problem; 
in practice these equations are solved with approximate methods, for example, by 
expanding the HFB wave functions in a harmonic oscillator basis [3] . 



68 



The Green's functions Go and Go are related to each other by Go(e) = 
-G* (~e*). 



Similarly, for G s and F we can write the set of two equations 

G s (e) = G (e) - G (e)A*F(e) , F(e) = G (e)AG s (e) . 



(B.5) 



From the Gor'kov equations it follows that 

G s (e) = -G s *(-e*) , F*(e) = -F*(-e*) . 



(B.6) 



The "zero" Green's function Go is diagonal in the basis of eigenfunctions of 
the Hamiltonian h (/^-representation), 



d2/i(l,2)0 fc (2) = e k (j) k (l), 



(B.7) 



and for its spectral decomposition we get 



G (l,2;e) = £ 



k l 



+ 



1 - n k 



e - e k + ji - i5 e - e k + jj, + i5 



J fc (l)&(2), (B.8) 



where n k = 1 if e k < jj, and n k = if e k > fi. 

For systems with even particle number, the spectral expansions of G s , F and 
F\ in coordinate space representation, may be written in the form: 



G s (l,2;e)= E 

(£ Q >0) L 



V:il)V a {2) | £f a (l)!7*(2) 
e + -E a — i<5 e — + i<5 



i ? (l,2;c)= E 

(£ a >0) 

i ?t (l,2;c)= E 

(£ Q >0) 



V:(l)U a (2) | ^(1)^(2) " 
e + — i<5 e — + i5 

^(l)K(2) | K(i)^(2) ' 

e + E a — i5 e — E a + i5 



(B.9) 
(B.10) 
(B.H) 



with E a the exact eigenvalues and U a , V a the exact eigenfunctions of eqs. 
(A. 28). The corresponding expressions in the /c-representation are: 



<?.y(e)= E 

(£ Q >0) L 



V* V- U- u* 



e + E a — iS e — E a + iS 



[B.12) 
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F^)= E 

(E a >0) 



V* U- 



+ 



U- V* 



e + E a — i5 e — E a + iS 



(B.13) 



4w= E 

(E a >0) I 



u ia v ja 
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V- u* 



e + E a — iS e — E a + id 



(B.14) 



with U ka = J dl 4>* k (l)U a (l) , and V ka = J dl fc (l)V Q (l) . In diagonal approx- 
imation, which is exact in uniform infinite matter, the only nonzero matrix 
elements of G s , F and F^ are: 



with I A;) being the time reversed of \k) as defined by (A. 62). Here we have 
used the customary conventions: 

A fc = -2E k v* k u k , E k = \/(e k -/i) 2 + |A fc | 2 , 
2 1 e fc - /i , |2 1 e k — /ix 

ki =2( 1+ ^r)> i^i = 2 (1 -^r } - 



The corresponding spectral expansions for G s may be found from those for 
G s by using the relations (B.6). It is easily seen that the expansions written 
above for the anomalous Green's functions F and F^ may be also obtained 
from each other by using (B.6). 

In what follows the Green's function method is used to extract the pairing 
contribution to the total energy of the system. To isolate the pure pairing 
part one has to consider two effects: (i) the direct influence of the pairing 
gap A at fixed single-particle Hamiltonian h — fi and (ii) the "polarization" 
mechanism due to variation of h and \x in the presence of the pairing field. 

To obtain the first change in the density, when A — > but h and /i fixed, we 
introduce a "zero" density matrix: 

p (l,2)= / ilGo(l,2;e). (B.18) 
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For the difference 5p = p — po, eq. (B.3) yields 

de „ „ „ x „ „ , s r de 







2 







Using (B.8) and (B.14) we obtain in the /c-representation: 

8 Pa = ~ E E F Ti N^^a + (1 - m)Utyia] , (B.20) 



or, in diagonal approximation, 



W r) = E E(! " H) oWP '^ _„J <M^ ^) I 2 • (B.21) 



One should keep in mind that p and po belong to states with different particle 
number: 

5N(p) = N(p)-N (p)^0. (B.22) 



The second change in the density we get in the situation without any pairing. 
Then we have the HF vacuum and the corresponding HF Green's function, 

^H F (e)= 1 , (B.23) 

e — AlHF + Phf 



with hu_F being the self-consistent HF single-particle Hamiltonian. The HF 
density matrix is given by 

P hf(1,2)= / ^1Ghf(1,2;6) (B.24) 



with the average number of particles (HF|7V|HF) = ./Vhf(phf), i-e. Trp H F = N. 

Collecting contributions from both pairing-induced effects we get the total 
variation of the density matrix: 

fypair = p - Phf = Sp + 5p , (B.25) 



with 5po = po — Phf- From the definition of the Green's function we find 

G HF - Go 1 = h - p - h HF + phf , (B.26) 
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or 



Co — Ghf — GW(<5t/ P air — S/I)Go , 

where SU pa i T = h — /ihf and 5p = p — /zhf- 
The expression for Spo can then be written as 

$Po = X^( G °(e)-G HF ( e )) 

= / ^G HF (e)(5U pair - 5p)G (e) , (B.28) 

Dealing with the first-order terms of perturbation theory, the integral over e 
may be replaced by the static particle-hole propagator A (see Ref.[10]): 

Ic ^[ G ^ G °^ ~ j c ^ G H F (e)G H F(e) =A. (B.29) 

Its convolution with any perturbation operator which is diagonal in the HF 
basis, i.e. commutes with /ihf, vanishes (no polarization effect). This property 
will be used to obtain the energy variation connected with 5p (see eqs. (B.38) 
and (B.39) below). The average number of particles of each kind is fixed: 
(HFB|iV|HFB) = (HF|JV|HF). The total density variation due to pairing, of 
course, does not change the particle number, so 

5iV pair = 5N(p) + 5N (p, p UF ) = . (B.30) 

By definition, the pairing energy is the total change of the system energy due 
to pairing correlations: 

-Kpair = ^normal [p] + -^anomal [p, v\ — -^normal [PHf] • (B.31) 

Analogously to the procedure applied above to the density matrix, we can 
split the pairing energy into two parts 

E pair = 5E + 5E , (B.32) 

where the difference 

8E = -EnormalM + E anoma l[p, v\ — -^normal [Po] 



(B.27) 
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is the "direct" pairing contribution and 

5Eq = -Enormal[po] — -^normal [PHf] 

is the polarization pairing energy. 

In the first step we calculate 8E. The anomalous part of the energy (A. 34) 
may be written in the form: 

£anomal[p, u] = ^Tl jf ^AF* (e) . (B.33) 

This can be expressed, by using the Gor'kov equations, as 

^anomai = ^Tr / ^(e - h + P )G s (e) . (B.34) 
2 Jc 2.H1 

On the other hand, the change of the normal part of the energy functional, in 
first-order perturbation theory, is given by 

^normal [p] ~ ^normal M = Tr(h5 P ) = fi5N + Tr[(h ~ fi)5 P ] . (B.35) 

This expression may be written in another form: 

^normalH - £„ormal[Po] = fl5N + Tr ^ ^(h - fi)(G s (e) - G (e)) . 

Adding -E an0 mah we find the first ("direct") change in the energy connected 
with the density variation 5p: 

5E = ptSN + ^Tr / ^(e + h- pt)G s (e) 

Z JC ZlTl 

-Tr / ^.(h-iM)G (€), (B.36) 
which can also be written, by using the Gor'kov equations (B.3), as 

SE = fiSN + ^Tr / ii( e + / l - Ai )G ( e )AFt(e). (B.37) 

2 Jc 2.TII 

For the second ("polarization") contribution to the total pairing energy in 
eq. (B.32), related to the density variation (B.28), we get: 
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SE = Tr(h KF 5p ) = Hb.fSN 

+ Tr / ^(/ihf - m F )G HF (e)(5U pair - 5^)G (e) . (B.38) 

The last term vanishes in first order. Thus, for the second change in the total 
energy, again in the first order, we find 

5E = fi HF 5N « fi5N . (B.39) 



Collecting both contributions and taking into account the particle conserva- 
tion (B.30), which gives fi(5N + SN ) = 0, we arrive at the resulting formula 
for the pairing energy: 

£ P air = ^Tr J^L( e + h- //)G (e)AFt(e) . (B.40) 



In /c-representation, using the spectral decompositions of eqs. (B.8) and (B.14), 
this can be written as 



E ■ 

-^pair 



1 y, yv E a I e j /i I y^ 



'J 



^ i (E a >0) Ea + l 6i ^1 j 

{(v ja u; a - u; a v ia ) - (1 - ^(v;-^ + c^vy 

In diagonal approximation we get: 



Ei 



pair 



-5E 



|A fe | 2 








2£ fc 









(B.41) 



^B.42) 



For infinite nuclear matter the sum over k is replaced by an integral in the 
momentum space, and one gets the pairing energy density (for one kind of 
particles): 

1 r dp | A (p) | 2 E(ft) - \e(ft) - e F 
£ P air = ~2 J ( 2vr )3 E (p) E(p) + \e(p)-e F 



(B.43) 



In the simplest model we assume 



|A(^)| = |A(pp)| = A, c(^) = |^» (B.44) 



which leads to the expression 



3 A 2 7 VTTt VWT¥-\t\ 
£ P air = —zQo — / dt — , (B.45) 

8 e F 7 V5 2 + t 2 V<5 2 + t 2 + t 



-1 
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with 

t = 1 , S = — , Q = —. 

The integral in (B.45) at 5 <C 1 can be easily estimated to yield 
1 — 5 2 (6 In 2 + 1 — 2 In 5) /32 + higher order corrections in 5. Numerically this 
integral deviates very little from unity even at unreasonably large 5 (less than 
by 2.5% up to 5 « 0.5). 

Thus, in the weak pairing approximation for an infinite system, as long as the 
condition A ep holds (which is generally assumed for nuclear matter near 
the saturation point), the pairing energy density is 

3 A 2 (p F ) lm*p F . . 

^pair = " 5 £0^^ = --^A 2 ( PF ) . (B.46) 
O €f 4 7T 

We want to conclude this Appendix with the remark that in the case of contact 
pairing force leading to a local pairing field A(r) there is no difficulty related 
with the cutoff e c to calculate the total energy of the system: in the energy 
space as a function of e c it converges rapidly (<5 C -E pa i r ~ l/tcy^) although 
kinetic and interaction energies, if taken separately, diverge (~ +^/^ and 
~ — y/e^, respectively). 



Appendix C. Solving the Gor'kov equations in the coordinate-space 
representation 

In this Appendix we describe the coordinate-space technique which allows, in 
the case of spherical symmetry, to solve the Gor'kov equations for given local 
mean-field and pairing potentials exactly. This technique has been invented 
in [25], where the most interesting effects, arising due to the proper treat- 
ment of the coupling between bound orbitals and particle continuum, were 
carefully investigated. Among them are the lowering of the chemical potential 
/i (an increase of the binding energy), the term-repulsion phenomenon and 
the appearance of width for deep-hole states lying below 2/x from the contin- 
uum threshold. It has then been applied, within the EDF approach with a 
density-independent pairing 5-force, to some superfluid even-even nuclei [20]. 

The fact that the coordinate-space HFB approach for finite systems can prop- 
erly treat the positive-energy continuum, in contrast to BCS-like methods in 
which the presence of an unphysical "particle gas" is almost unavoidable, was 
apparently first pointed out in [108]. It has been shown there that the HFB 
approach naturally leads to a localized nuclear wave function and gives the 
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correct asymptotics for the normal and anomalous densities. These properties 
of the HFB solutions are especially important for the nuclei close to the drip 
lines. A method of solving the HFB equation directly in the coordinate rep- 
resentation, but approximating the continuum by a set of the discrete states 
in a spherical box, has been used in [74] for the self-consistent description 
of a chain of tin isotopes with the Skyrme effective interaction. The asymp- 
totic behavior of the quasiparticle wave functions, as well as of the density 
distributions, have been studied with this method [74,109]. 

Here we give a more detailed description of the technique suggested in [25] for 
spherical even-even nuclei. We shall present also, by using the so-called "uni- 
form filling approximation" to preserve the spherical symmetry, an extension 
of this technique for the odd systems. With this extension the blocking effect 
appears in a natural way. In what follows, the isospin variables are omitted to 
simplify the notation, and the derivation is the same both for neutrons and 
protons. 

In order to separate the spin-angular variables in the Gor'kov matrix equation 
(B.l), it is convenient to introduce a unitary transformation for the generalized 
Green's function: 

G = f ] gf, g = fGf\ (c.i) 

where T is given by 

' 6(1, 2) 

T= \ K J . (C.2) 

T(l,2) 



The operator T(l,2) is defined by eq. (A. 63). Applying this transformation 
to (B.l), we get the equation for Q: 

(e-H)g = I, (C.3) 



where the transformed Hamiltonian reads 



n = fnf ] 



' h - n Art N 

v -fA* n-Th*t* j 

V,2) -^(1, 2) A(fi)*(l,2) 
v A*(fi)*(l,2) ^(1,2) -Ml, 2) 



(C.4) 
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Here time-reversal invariance, Th*T^ = h, is assumed (no magnetic field or 
rotation) and the local (real) pairing field A(r), which is diagonal in spin 
variables, is introduced according to eq. (A. 66). The transformed generalized 
Green's function Q entering eq. (C.3) is given by 



0(1, 2; e) 




(C.5) 



Because of the time-reversal symmetry, required for the Bogolyubov quasipar- 
ticle vacuum, this transformed Green's function may be written, for even-even 
spherical nuclei, in a form where the spin-angular parts are the same for all 
its 2 x 2 components: 

2 ; e ) = — Yl 9ji( r u r 2] e)0j7m(ni, si)0*, m (n 2 , s 2 ) , (C.6) 

r l r 2 jlm 



with (f)ji m (n, s) = Cj™_ s i i Yj m _ s (n) being the s-component of the usual spher- 
ical spinors (s = ±|). For the generalized radial Green's function gji, which 
has the matrix form 



9ji(r,r';e) 



9ji(r,r';e) g}f(r,r';e) 
gf(r,r';e) g]f(r,r';e) 



(C.7) 



one gets the equation 



e - hji + ji -A 



v -A e + hji-n^ 



g j i(r 1 ,r 2 ;e) 



6(n - r 2 ) 

S(n - r 2 ) J 



,(Ci 



h 2 / d 2 , 1(1+1)- 



+ Uji(r) is the single-quasiparticle Hamiltonian 



where h jX = ^(-£2 
in the jl channel (see below). 



The solution of the matrix equation (C.8) can be constructed by using the set 
of the four linearly independent solutions 



ViA r ) 



u i,jl( r ) 

Vi,ji(r) 



, 2 = 1-4, 



which satisfy the homogeneous system of equations obtained from (C.8) by 
setting the right hand side to zero: 
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( h 2 d 2 \ 



A ( r ) u i,ji( r ) 



2m dr 2 



fji(r)+ti-e)v iijl (r) = 0. 



(C.9) 



Here 



f jl (r) = U c (r)+U al {r)(ffl) jl + 



h 2 i(i + i) 

2mr 2 



U c and U s i are the central and spin-orbit mean-field potentials, respectively, 
= j(j 1 + 1) — 1(1 + 1) — f , m is the nucleon mass. We shall seek the radial 
Green's function gji entering eq. (C.6) in the following form (for simplicity, 
the e-dependence and the jl indices are not shown): 



2m 



g(ri,r 2 ) = -p-[(Q!J/i(ri) + PhiXi)) y 3 (r 2 ) 

+ (iyi( r i) + ^2(^1)) y^r 2 )]0(r 2 - n) 

+ -^[fc^i) (02/1(7-2) + Py 2 (r 2 )) 

+ yi(n) (7^1(^2) + Sy 2 (r 2 ))}6(ri - r 2 ) , 



(CIO) 



where y i = (ut ,Vi) is the transposed i-th solution, 9(x) is the step function, 
a, P, 7, and 5 are the coefficients to be determined. The solutions y\ and y 2 
are chosen to be regular at r — > and the other two, y 3 and j/ 4 , to be regular 
at r — > 00: 



yi(r-»0) 



y 3 (r -kx>) 




( g+r ) ,+1 , t/ 2 (r^0) 



2/4 (r ^00) 




(g_r) 



(C.ll) 



Here 



5± 



2;/) 



A 2 
1/2 



nl/2 



A 2 
^0 



Co = -A / ^ + v /e 2 

(/xi^A 2 ^]^, Coo = -A oc /(e+ V^A 2 :) ; ( C - 12 ) 
U c (0) + ^(0)( c ?r) , A = A(0) , Aoo = A(r^oo) . 



It is important to notice that on the physical /c-sheet the imaginary parts of 
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the asymptotic momenta k± should be chosen positive^]: Imk± > 0. 



Using the system of equations (C.9), it can be easily verified that for any two 
arbitrary solutions, 



and jjj 



the combination Cjj = U%j — V%j composed of the Wronskians U, 
and Vij = ViVj — v^Vj, where the primes denote differentiation with respect 
to r, does not depend on the coordinate r, and for solutions satisfying the 
boundary conditions (C.ll) one gets Cu = C34 = 0. It can be seen also that 
U13U24 — U23U14 = U12U34 and V13V24 — V23V14 = V12V34. With these relations 
one easily obtains an equation for determining the four coefficients entering 
eq. (CIO): 






(C.13) 



where D(e) = C13C24 — C23C14 is the determinant of the matrix C . 

The formulas (C.l), (C.6) and (C.10)-(C.13) constitute an exact solution of 
the Gor'kov equation for the generalized Green's function in the coordinate 
representation for superfluid spherical nuclei with a localized pairing field 
A(f). 

In the case of vanishing pairing one has w 2 = w 4 = v\ = V3 = and (3 = 7 = 0. 
Then the radial generalized Green's function (C.7) becomes diagonal with g 11 
and g 22 in the following form: 



.11 A_ 



2mu 1 . j i(r < ;e)u 3;j i(r > ;e) 



J h U 13 {e) 

22, „/. _ 2m v 2;jl (r < ;e)v 4 ., j i(r > ;e) 

9jl {r,r,e)\ A=0 - ^ , (C.14) 

where r< and r > denote the lesser and the greater of r and r', respectively. 
Here g 11 is interpreted as the HF Green's function for particles (e > fi) while 



16 To be more precise, the solutions ^3 and 2/4 are given asymptotically at r — > 00 by 
the Whittaker functions W / _ i ^ ± ^ +1 / 2 (— 2ik±r), respectively, with rj± = mZe 2 /h 2 k± 
being the corresponding Coulomb parameters {rj± = for neutrons). 
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g 22 as that for holes (e < /i)^ 7 ]. By integrating the latter along the contour 
C of Fig. 21 one gets the HF jl density matrix Pji(r, r')|A=o which is normal- 
ized by (2j + 1)/ drpji(r,r')\/^ =0 = Nji where Nji is the number of particles 
occupying orbitals with the given quantum numbers jl. This density matrix 
can be expressed by a finite sum of separable terms, i.e. by a sum of residues 
of <7 22 |a=o at the poles which are zeros of the Wronskian V24. The equation 
V24 = defines the spectrum of the HF hole states, and this spectrum is 
utterly discrete. 



Im 



-<-F 



Re e 



Fig. 21. The contour C in the energy plane to integrate the generalized Green's 
function for determining the generalized density matrix (for even nuclei) . The crosses 
mark the positions of the quasiparticle poles, the heavy lines represent the branch 
cuts on the Ime = axis where the energy spectrum is continuous, fi is the chemical 
potential, €f is the Fermi energy and e c is the energy cutoff (see text). 

The situation with pairing is completely different. It is clear that the poles 
(if any) of the generalized Green's function are determined by the equation 
D(e) =0 which gives the discrete spectrum of the Bogolyubov quasiparticle 
states. Due to the invariance of the equation (C.9) under simultaneous re- 
placements e — > — e and (u, v) — > (v, —u) the energy spectrum is symmetrical 
around the point e = 0. For bound systems, by definition, one has \i < 0. 
If the system is finite, i.e. if the density distribution has a finite range, one 
expects Aqo = 0, otherwise the system would be unstable with respect to two 
particle emission. From eqs. (CIO)— (C.13) it follows then that the quasipar- 
ticle spectrum is discrete within the energy region |e| < (here gji is a real 
function at the axis Ime = 0) and continuous if |e| > (there §ji is a com- 
plex function). These features are reflected in Fig. 21 where the branch cuts 
are shown by the heavy lines extending symmetrically to the left and to the 
right from the points ±/i, respectively, along the Ime = axis. The case with 
only the branch cuts, without poles, corresponds to a drip-line even nucleus. 
Performing the integration of the radial Green's function §ji along the contour 



17 Eq.(C.14) for the single-particle Green's function allows one to construct the 
particle-hole propagator in the coordinate representation and to solve, for closed- 
shell nuclei, the RPA equations exactly, including the whole particle continuum. An 
earlier application of this method for the study of pion condensation in finite nuclear 
systems may be found in [110], and, for calculating the continuum-RPA multipole 
response functions, in [111]. 



80 



C of Fig. 21 one obtains the radial part of the generalized density matrix^]. 
To calculate the normal and anomalous densities we need the following two 
radial u jl density matrices": 

Pji(r, r') = J c ^f(r, r'; e) , ^(r, r') = jf ^^ 2 (r, r'; e) . (C.15) 

Then the resulting expression for the normal density reads 

P(f) = iD 2 J + l)P*i(r,r) • (C.16) 



It is normalized by / dfp(r) = A" which is actually the condition (2.23) for the 
chemical potential /i. For the anomalous density from (A. 37), (A. 64), (C.5), 
(C.6) and (C.15) we get 

K*=) = xLE0" + i)^(r,r). (C.17) 



These densities are used to calculate both the normal and anomalous parts 
of the interaction energy in eq. (2.22). In addition, for the spin-orbital term, 
eq. (3.5), we define^} 

= T^E( 2 i + l)(^~W(r,r) . (C.18) 
7rr jl 

18 The analytical properties of the Green's function allow one to consider only the 
upper half of the contour with Ime > 0: integrate ReG along the vertical sections 
and ImG along the horizontal section in the upper half-plane, then for the total 
result take the doubled sum of these integrals. 

19 Remember that we do not show here the cutoff index c. It should be emphasized 
that because of the continuum all the densities defined in this Appendix cannot 
be expressed by a finite sum of separable terms. The latter would be possible only 
within a certain approximation, namely with a discretized basis and with an energy 
cutoff. The coordinate-space technique described above, which involves integration 
of the exactly constructed generalized Green's function in the complex energy plane, 
appears to be exact for the j7-densities: both Pji(r,r) and Vji(r,r) converge as 

functions of e c (5 c Pji ~ e^ 1 an d S c i/ji ~ e c 1//4 , respectively). It is exact also for 

— 1/2 

the normal density C.16) since it converges, 5 c p ~ e c . However, if one goes 
beyond the cutoff e c , the anomalous density (C.17) diverges, v ~ e^ 4 . The kinetic 
and anomalous energies beyond e c diverge too, but their sum converges very rapidly 
((Jc-Epair ~ e c , see Appendix B). This fact has been used in the present paper to 
formulate the generalized variational principle with the cutoff local EDF. We stress 
once more that this formulation does not imply a cutoff of the basis. 
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The kinetic energy of the system is given by 



d 2 1(1 + 1) ' 



dr 2 



Pji\T,r 



(C.19) 



Note that to compute this energy one needs not only the normal ^/-densities, 
but also the jl normal density matrices pji(r,r'). 

For an odd system the spectral expansion of the generalized Green's function, 
eqs. (B.9)-(B.10), should be modified. Suppose that the addition of an odd 
particle leads to the appearance of the quasiparticle with the energy E aQ in the 
ground state of the system. Here cto denotes all the relevant quantum numbers. 
They can be specified, for example, by «o = {njlm}o from the whole set {a} 
used to label the solutions of the HFB equation (A. 28). Suppose further that 
E a corresponds to a certain eigenenergy of this equation and belongs to the 
discrete spectrum, i.e. > E ao > 0. As illustrated in Fig. 22, E ao is located 
on the Ime = axis in the vicinity of the point e = 0. Generally, E ao is of 
the order of A, the average matrix element of the pairing potential on the the 
Fermi surface. Note that the pole nearest to e = need not correspond to the 
ground state of the system. The case of E ao ps —/i determines the position of 
the drip line for odd nuclei. The rules for passing the Green's function poles 
for superfluid odd systems are formulated by Migdal [10]. To the spectral 
expansions (B.9)-(B.10), which do not contain the contribution of the odd 
quasiparticle, one should add the following terms: 



5Gf d (l,2;e) 



V*(l)V ao (2) , U ao (l)U* (2) 



"0 

e + E ao - ib 



+ 



Q[) ' 

e - E ao - ib 



| ^(1)^(2) | ^(1)^(2) 
e + Ea n + id t- E^, + ib 



(C.20) 



V*(l)U ao (2) , U ao (l)V*(2) 

Or l i ,2;ej = — — — H — 

e + E ao - id e - E ao - lb 

vim^m UaoWzm (C2l) 

e + E ao + ib e-E^+ib' 

The subscript a in these expressions, in Migdal's terminology, refers to an 
odd quasiparticle added to a system of N particles and to an odd quasihole in 
a system of N + 2 particles. Consequently, the subscript d refers to the odd 
quasihole in the system of N particles and to the odd quasiparticle added to a 
system of iV — 2 particles. Thus the four terms in (C.20) and (C.21) correspond 
to pairs of time-reversed states. But, due to the signs of the infinitesimal 
imaginary parts +ib in the denominators, only two, not time-reversed, terms in 
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Im e 



Q 



a 



Re e 



Fig. 22. The contour C for odd nuclei. The positions of the poles for an odd quasi- 
particle with the energies iE ao are shown by heavy dots. Note that the ways of 
passing these poles go in opposite directions. Other notations are the same as in 
Fig. 21. 

these expressions should contribute to the generalized density matrix. Closing 
the contour in the upper half of the energy plane, the terms from the first lines 
in eqs. (C.20) and (C.21) will be accounted for. By integrating the Green's 
functions G s and F along the contour C of Fig. 21, all the residues at the 
pole e = — E ao (= —Ea ) have been already included. Therefore, from the so 
obtained density matrix we have to subtract the contributions related to the 
third terms of these equations and to add the residues from the second ones. 
The resulting contour for the odd system is depicted in Fig. 22 which clearly 
illustrates the blocking effect: the presence of the odd particle in the level «o 
prevents it from participating in the pairing correlations because in this case 
the "conjugate" level d should be empty. This way we get the contributions 
to the normal and abnormal density matrices from the odd qusiparticle: 

5p odd (l,2) = -^ o (l)\4 (2) + c/ ao (l)f/; ) (2), (C.22) 

^° dd (l,2) = -^ o (l)^ () (2) + c/ ao (l)K: o (2). (C.23) 



In the following we shall consider the case when only the spherical part of 
the core polarization induced by the odd quasiparticle is important, and the 
deformation of the core can be neglected. Then the solutions of the Bogolyubov 
equation may be chosen to have the form 

U a (f, s) = ^u nj i(r)(j) jlm (n, s) , 

K(r, s) = ^„>)(-ir m+ V^_ m (n, s) , (C.24) 

where the radial wave functions u n ji(r) and v n ji(r) satisfy eq. (C.9). For the 
discrete states, including the odd quasiparticle state with a = ao, the latter 
obey the orthogonality relation / dr[u*, y (r)u n y(r)+?;*,y(r)?; n y(r)] = S nn >. We 
introduce further the uniform filling approximation. Namely, we assume that 
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the odd quasiparticle occupies all the possible tjiq states with equal probability 
and take the averages 



5p odd (l,2) = — --5> 0dd (i>2), 

6^(1, 2) = -J— £ ^ odd (l, 2) . (C.25) 

470 + J- ra „ 



Such an averaging procedure, with the wave functions given by (C.24), enables 
us to separate the spin-angular variables for the odd spherical system in the 
same way as done above for even systems. The resulting formulas for the 
contribution to the radial jl density matrices (C.15) due to the presence of 
the odd particle read 



< 5 /CL( r ' r ') = -<ojo«o( r )*Wo( r ') + u n jolo( r ) u *nojolo( r ') ' 

<5>, r') = -v* nojolo (r)u nojolo (/) - u nojolo (r)v* nojolo (r f ) . (C.26) 

The contributions to the normal and anomalous densities (C.16) and (C.17) 
are given by 

^5>1 = -Ar 2 {-Koioio(r)\ 2 + \u nojo io(r)\ 2 ) ■ (C27) 



Attt 2 

4nr 2 



5 <lio (0 = xL^Wo (r)v* nojolo (r) . (C.28) 



It can be seen that 8v° d j olo (r) has the opposite sign as u{r), eq. (C.17), reflect- 
ing the fact that the blocked level cannot directly contribute to the anomalous 
energy. The nucleon separation energies S n are determined by the position of 
the poles close to e = 0. As easily understood from Fig. 22, since these separa- 
tion energies are measured from the continuum threshold e = — Li, for an odd 
system one gets S° dd —li — E ao , and, for an even system, S^ cn —li + Ea . 
Thus we have S° ven - S° dd « 2E ao w 2A, i.e. the familiar odd-even effect in 
nuclear masses. 

The rest of this Appendix is devoted to the discussion of the asymptotic be- 
havior of the densities p{r) and u{f). Consider first the asymptotic properties 
of the wave functions y^ used to construct the radial Green's functions (C.10) 
through which these density can be obtained. For finite systems, the functions 
v%{r) and w 4 (r) in the solutions y^ and y± of the radial HFB equation (C.9), 
as given by the boundary conditions (C.ll), vanish at large distances since 

Coo ~ ^oo 

= 0. The function « 3 at large r behaves as exp ik + r and either 
decreases ~ exp(—r\j2m(\e\ + fx)/h 2 ) if |e| < —li or oscillates if |e| > —fx, 
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the latter corresponds to the continuum spectrum. The function V4 behaves 
asymptotically as expi/c_r. If the system is bound, i.e. fj, < 0, this function at 
r — > 00 always decreases exponentially ~ exp(— r^2m{\e\ — p)/Ti 2 ). 

To get the asymptotic laws for normal and anomalous densities, defined by 
eqs. (C.15)-(C.17), we have to consider now the radial Green's functions 
g 12 oc (7M1 + 5u2)v4 and g 22 oc (71*1 + <5t>2)t>4. Asymptotically, at large r, 
the functions Uj and Vi, i = 1,2, which are regular at r — > 0, are given by 
Ui(r) = q exp(i/c+r) + di exp(— ik+r) and i>i(r) = ejexp(i£;_r) + /, exp(— i/c_r), 
respectively. Using eq. (C.13) one can easily verify that, by construction, the 
combination 7^1 + 5o?2 vanishes. The combination 7/1 + 5/2 vanishes at the 
poles of the Green's function determined by the equation D(e) = 0. Finally, 
as the result of of integration of the Green's function along the contour C, we 
get 



p(r -> 00) ~ exp(-2r^/2m(|e | - /i)/ft 2 ) , 

v{r -> 00) ~ ^ exp(-r[ v /2m(-|e | -/i)A 2 + ^2m(|e | - ^)/^ 2 ]) , (C.29) 

where eo is the position, within the contour C, of a singularity point nearest 
to e = 0. This point may be a pole of the Green's function or the beginning 
of the branch cut. The latter situation corresponds to the drip-line nuclei. 
Generally, |eo| is of the order of the average matrix element A of the pairing 
potential on the Fermi surface (typically, A 1 MeV). In stable nuclei, since 
the absolute value of the chemical potential is much larger, // —8 MeV, 
the asymptotic behavior of both densities (C.29) differ but insignificantly. 
Approaching the drip-line, in even system, |eo|, and A may become of the 
same order: |e | ~ —fJ> ~ A. In such a case one gets 



p{r — > 00) ~ — exp(— Ar\J mA/h 2 ) , 



r 



v{r -> 00) ~ — exp(-2r v /mA/r) . (C.30) 



It follows that in this situation the anomalous density has a longer tail with 
decaying length by a factor of two greater than that of the normal density. 
Similar conclusions have been drawn in [74,108,109]. On the other hand, in 
the genuine drip-line even nuclei, with only one bound state - the ground 
state, when all the HFB eigenstates are embeded in the continuum, there are 
always contributions to the densities from the quasiparticles with different 
quantum numbers jl and with the energies close to /i. The weights of these 
contributions are determined by the corresponding quasiparticle level density, 
but in any case the most distant tails of the densities in drip-line nuclei are 
asymptotically given by the expressions 
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p(r -> oo) ~ — exp(-Arym\fi\/h 2 ) , 

v{r -> oo) ~ exp(-2ry / m|/i|/^ 2 ) , (C.31) 

which are obtained from eq. (C.29) by taking the limit | e 1 ~~ *■ \iA- The tech- 
nique described in this appendix permits, in principle, an exact calculation of 
the asymptotic behavior of the HFB densities in different situations. 
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